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SUMMARY 

Panel methods are numerical schemes for solving (the 
Prandtl-Glauert equation) for linear, inviscid, irrotational 
flow about aircraft flying at subsonic or supersonic 
speeds. The tools at the panel-method user’s disposal are 

(1) surface panels of source -doublet-vorticity distributions 
that can represent nearly arbitrary geometry, and 

(2) extremely versatile boundary condition capabilities 
that can frequently be used for creative modeling. This 
report discusses panel-method capabilities and limitations, 
basic concepts common to all panel-method codes, differ- 
ent choices that have been made in the implementation of 
these concepts into working computer programs, and var- 
ious modeling techniques involving boundary conditions, 
jump properties, and trailing wakes. An approach for 
extending the method to nonlinear transonic flow is also 
presented. 

Three appendixes supplement the main text. In 
appendix A, additional detail is provided on how the basic 
concepts are implemented into a specific computer pro- 
gram (PAN AIR). In appendix B, we show how to evalu- 
ate analytically the fundamental surface integral that 
arises in the expressions for influence-coefficients, and 
evaluate its jump property. In appendix C a simple exam- 
ple is used to illustrate the so-called finite part of 
improper integrals. 

1. INTRODUCTION 

Panel methods are numerical schemes for solving (the 
Prandtl-Glauert equation) for linear, inviscid, irrotational 
flow at subsonic or supersonic free-stream Mach numbers. 
Currently, panel-method codes are the only codes com- 
monly in use that are sufficiently developed for routinely 
analyzing the complex geometries of realistic aircraft The 
objective of this report is to give the reader some idea of 
what panel methods can and cannot do, to describe their 
common roots, to describe the differences between 
various specific implementations, and to show some 
example applications. In addition, recent progress in 
solving nonlinear transonic flow problems by combining 
portions of panel-method technology with other numerical 
techniques, is described. This material is followed by 
three appendixes that contain additional details: appendix 
A describes how the basic ideas common to all panel 
methods are actually implemented in a specific code; 
appendix B shows how to evaluate some of the integrals 
that arise in the influence-coefficient equations; and 
appendix C discusses the so-called finite part of improper 
integrals. 


There are fundamental analytic solutions to the 
Prandtl-Glauert equation known as source, doublet, and 
vorticity singularities. Panel methods are based on the 
principle of superimposing surface distributions of these 
singularities over small quadrilateral portions, called 
panels, of the aircraft surface, or to some approximation 
to the aircraft surface. The resulting distribution of super- 
imposed singularities automatically satisfies the Prandtl- 
Glauert equation. To make the solution correspond to the 
desired geometry, boundary conditions are imposed at 
discrete points of the panels. (Mathematicians refer to 
these discrete points as collocation points; panel-method 
users refer to them as control points.) 

Panel codes are often described as being lower-order 
or higher-order. The term lower-order refers to the use of 
constant-strength singularity distributions over each panel, 
and the panels are usually flat. Higher-order codes use 
something of higher order than constant, for example, a 
linear or quadratic singularity distribution, and sometimes 
curved panels. 

Panel methods were initially developed as lower- 
order methods for incompressible and subsonic flows 
(e.g., refs. 1, 2; see ref. 3 for a review of panel methods 
existing through about 1976). The first successful panel 
method for supersonic flow became available in the mid- 
1960s (refs. 4, 5). This was also a lower-order method, 
and is variously referred to as the constant-pressure panel 
method, or the Woodward-Carmichael method. 

The panel methods for three-dimensional subsonic 
flow allowed the actual vehicle surface to be paneled, 
whereas the Woodward-Carmichael method was more 
severely restricted in the placement of the panels. For 
example, the wing was a planar array of panels, the body 
(fuselage) volume was modeled with a line distribution of 
source and doublet singularities (resulting in a body-of- 
revolution) and the body boundary conditions were 
imposed with a cylindrical “interference” shell of wing- 
type panels. This representation was later extended to 
include multiple wing-body components (ref. 6), but was 
still restricted to the planar panel representation. These 
two extremes of actual-surface models and mean-surface 
models (Woodward-Carmichael) are illustrated in fig- 
ure 1. 

The mean-surface model used in the Woodward- 
Carmichael panel method was a consequence of numerical 
stability problems that arose in supersonic flow. The con- 
stant-strength, elementary horseshoe vortex singularity 
distribution (producing a constant pressure over each 
panel) often produced unstable numerical behavior (the 
solutions would “blow-up”) when a panel was inclined to 



a supersonic flow. The method worked only when all the 
panels were kept parallel to the free-stream flow. This 
required that angle of attack, wing thickness, camber, and 
twist be simulated through the boundary conditions; that 
is, it was necessary to have the panels generate flow per- 
pendicular to themselves and thereby turn the flow 
through the desired angles, as is (tone in classic thin-wing 
theory. As a consequence of this restricted geometric 
model, several new approaches to the supersonic problem 
were pursued in the 1970s. 

The first of these was also due to Woodward; it 
evolved into the series of computer programs known as 
USSAERO (ref. 7). For fuselage panels, USSAERO uses 
constant-strength source singularities. Wing panels use 
elementary horseshoe-vortex singularities whose strength 
distribution varies linearly in the chordwise direction and 
is constant in the spanwise direction. Although this repre- 
sentation gave an improved modeling capability, numeri- 
cal problems would still often occur when the wing panels 
were inclined to a supersonic free-stream flow. 

Another approach, developed by Morino and his 
associates, uses a superposition of constant-strength 
sources and doublets on hyperboloidal panels. The con- 
stant strength doublets produce a velocity field that is 
identical to that produced by line-vortex elements, having 
the same strength as the doublet panel, placed head to tail 
around the panel perimeter (a so-called ring- vortex panel). 
This method is available in the computer program called 
SOUSSA; it too is unable to handle the steady supersonic 
case (ref. 8, pp. 2-20; private communication, L. Morino, 
Feb. 1981). 

The key to eliminating the numerical stability prob- 
lems associated with supersonic flow, was to use doublet 
distributions that were continuous over the entire surface 
of the aircraft. This approach, using quadratic doublet 
distributions (equivalent to linear vorticity distributions) 
was first used in the PAN AIR code (refs. 9-14) and its 
pilot code predecessor (ref. 15). It has since been imple- 
mented in the European version of PAN AIR, called 
HISSS (ref. 16). The continuous-doublet distribution 
eliminates the appearance of spurious line-vortex terms at 
the panel edges, which was the cause of the numerical 
stability problems in the earlier approaches. 

Within the limitations of the Prandtl-Glauert equa- 
tion, the higher-order singularity distributions used in the 
PAN AIR and HISSS codes allowed the actual-surface 
paneling models, long in use for subsonic flow, to also be 
used for supersonic flow. It also had a very beneficial side 
effect: the numerical solutions turned out to be much less 
sensitive to the size, shape, and arrangement of the panel- 


ing than in earlier methods, including the subsonic-only 
methods. Partly for this reason, continuous quadratic dou- 
blets were incorporated into the subsonic-only MCAERO 
code (ref. 17). These advantages did not come free how- 
ever. The higher-order distributions require much more 
analytic work to derive the influence-coefficient equa- 
tions, and demand many more arithmetic operations than 
the simpler lower-order (constant-strength) methods, 
which results in significantly higher run costs. 

It was subsequently discovered that for subsonic 
flow, setting the perturbation potential to zero at the inte- 
rior side of panels, in conjunction with the original lower- 
order singularity distributions, also reduced the solution 
sensitivity to variations in panel layout. This led to a 
renewed interest in the lower-order methods, resulting in 
the VSAERO (refs. 18, 19) and QUADPAN (refs. 20, 21) 
codes. QUADPAN was later revised to handle the super- 
sonic case by changing its constant-strength doublets to 
continuous linear doublets. 

2 . WHAT PANEL METHODS CAN AND 
CANNOT DO 

Panel-method-based computer programs are currently 
the workhorse codes for predicting the aerodynamics of 
complete configurations. Representative aircraft examples 
that have been analyzed with panel method codes are 
shown in figure 2. Although such codes are routinely used 
to analyze very complicated geometries, they do so at the 
expense of ignoring much fluid physics. The equation that 
panel codes solve is the Prandtl-Glauert equation. For 
steady subsonic flow this equation is usually written as 

V 2 <t» = (l-Mi)4» XA + «t)y y + <(>^=0 (1) 

and for supersonic flow it is sometimes multiplied by -1, 
-V 2 <t> = (Mi - l)<|> xx - $yy - ♦„ = 0 (2) 

where is the free-stream Mach number and <t> is the 
perturbation velocity potential. 

For subsonic free-stream flow, equations (1) and (2) 
are elliptic, being similar to Laplace's equation. Such 
equation types have the property that any disturbance at 
some point is felt everywhere in the flow field (although 
the effect usually dies out rapidly with distance). For 
supersonic free-stream flow the equations are hyperbolic, 
with the x-derivative term behaving like time in the wave 
equation. Solutions for the supersonic case are 
fundamentally different, disturbances having restricted 
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zones of influence (or in Von Karman’s words, zones of 
“silence,” or “forbidden signals”; ref. 22). The dis- 
turbances only propagate downstream, along rays defined 
by the Mach cones (characteristic surfaces), reflecting off 
downstream geometry and interfering in a wave-like 
manner with other disturbances. 

The Prandtl-Glauert equation is the simplest form of 
the fluid-flow equations that contain compressibility 
effects (i.e., the effect of Mach number on fluid density). 
It is obtained from the more general Navier-Stokes equa- 
tion by (1) neglecting all the viscous and heat-transfer 
terms; (2) assuming that the flow is irrotational, thereby 
admitting the introduction of a velocity potential; and 
(3) discarding all nonlinear terms. This restricts the flow 
to be inviscid, irrotational, and linear. Often, the flow is 
also assumed to be steady. Physically, these restrictions 
mean that important flow behavior such as separation, 
skin-friction drag, and transonic shocks are not predicted 
with panel methods. Items that are predicted include drag- 
due-to-lift (often called induced drag for subsonic flow, 
and vortex drag for supersonic flow), and wave drag. 

Wave drag is predicted because the Prandtl-Glauert 
equation admits solutions that approximate the weak- 
shock solutions of shock-expansion theory (ref. 23, 
pp. 215, 216). A simple example is the supersonic flow 
over a thin wedge (fig. 3(b)). For small wedge (deflection) 
angles, the shock is attached at the wedge leading edge, 
forms at an angle very nearly to that of the Mach angle, 
and the flow remains supersonic on the downstream side 
of the shock. The limiting case for these weak shocks, in 
which the shocks form at exactly the Mach angle, is pre- 
dicted by the Prandtl-Glauert equation. 

The absence of any explicit viscous effects causes 
subsonic flow solutions to be non-unique unless a Kutta 
condition at sharp trailing edges is somehow imposed 
(ref. 24, pp. 80, 81). This is done with the addition of 
some type of wake panels that trail downstream from 
lifting-surface trailing edges (fig. 3(a)), causing the flow 
to separate smoothly from these edges and allowing the 
potential to jump (be discontinuous) across the wake. 
Most panel methods require the user to assume the shape 
and position of the wakes. For a simple wing body this 
poses no difficulty, the wake position being relatively 
unimportant. However, for multiple-lifting-surface config- 
urations, the wake placement is important since it affects 
the flow experienced by downstream geometry. A few 
codes iteratively solve for the wake shape and location. 

Because panel methods are able to treat complete 
configurations, they have often been used in combination 
with other methods to approximately account for addi- 


tional physics neglected by the Prandtl-Glauert equation. 
One fairly common practice is to include the presence of 
the wing boundary layer (ref. 25). The basic idea is to use 
the pressure distribution from the panel -code solution as 
input to a boundary-layer code and compute the displace- 
ment thickness. This incremental thickness is then repre- 
sented in a second run of the panel code. This is usually 
done in one of two ways, as illustrated in figure 4 
(ref. 25). The first is to actually recompute the wing 
surface coordinates and the new wing-body intersection 
by adding the displacement thickness to the actual wing 
geometry. An alternative approach is to use “blowing,” in 
which the source strengths of the wing panels are adjusted 
such that each panel ejects (or sucks) enough fluid to 
cause the resultant flow field to be approximately 
displaced by the displacement thickness. For either 
approach, the resulting change in actual or apparent wing 
shape has two effects: it reduces the effective camber of a 
cambered wing and it increases the wing thickness. For a 
specified angle of attack, the primary aerodynamic effect 
of these changes is a reduced lift owing to the reduced 
camber. The second, but usually less important effect, is a 
slight increase in lift owing to the increased wing 
thickness. 

Another example is the coupling of panel codes to 
propulsion codes. In reference 26, the PAN AIR code is 
coupled to a parabolized Navier-Stokes propulsion code. 
The purpose was to account for the viscous, high-energy, 
exhaust-flow effect on the aerodynamic flow about the 
complete aircraft. 

Panel-method codes have also been built to model the 
flow separation that occurs off highly swept wings with 
sharp leading edges (refs. 27, 28). In these codes, wake 
panels emanate from the wing leading edge, as well as 
from the trailing edge (fig. 5). Iterative techniques are 
used to solve for the correct shape and position of the 
leading-edge wake panels. The criteria to be satisfied are 
(1) that the Kutta condition be enforced and (2) that the 
entire wake surface be a stream surface (i.e., no flow 
crosses it, and it supports no pressure jump). 

3. COMMON ROOTS OF PANEL METHODS 

As indicated in section 2, panel methods rely on sur- 
face distributions of sources, doublets, and vorticity. We 
will see later that doublet and vortex distributions are 
related. Since surface vorticity is a vector and a doublet is 
a scalar, it is often easier to work with doublets than with 
vorticity, and then compute the vorticity from the doublet- 
strength distribution. Most higher-order panel-method 
codes take this approach. 


3 



It can be verified by direct substitution, that the fol- 
lowing expressions, called unit point sources and dou- 
blets, respectively, satisfy the Prandtl-Glauert equation 
(eqs. (1) or (2)). 


This form clearly shows the directional properties of a 
point doublet and reveals that a doublet disturbance dies 
off at least as rapidly as the inverse of the distance 
squared. 


Point source: 



( 3 ) 


Point doublet: 


♦?(SqM-V q 


1 

r(x p ,x q ) 



( 4 ) 


where the so-called hyperbolic distance R is given by 
r = -\]( x q- x p) + P 2 (yQ-yp) + ( z q _z p) (5a) 


where 


P 2 =i-Mi 


(5b) 


In these expressions, point P is the influenced point in 
space having coordinates x P = (x P ,y P ,z P ) and point Q is 
the influencing point Xq = (xq^q.zq) at’ which the unit 
point source or doublet is located (see fig. 6). There is an 
elemental area dSQ associated with the doublet, and the 
doublet axis is normal to this area. (Recall from elemen- 
tary fluid mechanics that a doublet can be thought of as a 
source-sink pair approaching each other along an axis. 
This definition of a doublet produces the same result as 
eq. (4),) The subscript Q on the scaled gradient operator 
means that the derivatives are to be taken with respect to 
the coordinates of point Q, not point P. 

For incompressible flow, R becomes simply the geo- 
metric distance between the two points P and Q. Equation 
(3) then tells us that a point source at Q produces a distur- 
bance at P that diminishes inversely as the distance 
between the two points. The meaning of equation (4) is 
not so obvious until one works out the expression indi- 
cated by the dot product If one chooses the xyz coordi- 
nate system at point Q as shown in figure 6, then the unit 
normal n equals the unit vector k and equation (4) 
becomes simply 



-sinO 


( 6 ) 


Since the Prandtl-Glauert equation is a linear partial 
differential equation, sums of the source and doublet 
solutions are also solutions. Thus, panel methods are usu- 
ally thought of as superposition methods, and, hence, are 
restricted to linear problems. There is a more general 
approach, however, that, while containing superposition 
as a special case, can also be used to solve nonlinear prob- 
lems. In section 6 we will take a look at how panel- 
method technology can be combined with other tech- 
niques to solve the nonlinear full-potential equation, so it 
is advantageous to look at this more general approach, 
known as Green’s third theorem (ref. 29, p. 21 , eq. (7)). 

In reference 29, the derivation corresponds to incom- 
pressible potential flow; in reference 30, this is general- 
ized to the compressible case. The result is the following 
identity: 

<t>p - JJ^[<^(^q) k o(^q^p) + h(xq)k^(xq,x p )|iSq 

+JJJ v (V 2 <t>) K o(^x P ) dv (7 > 

where 

k o(*q.*p) = y<1>p(xq) 

k ( x(xq.xp) = 7«I>p(xq) 

K a (x,Xp) = -£-<t>p(x) 
dV = dx dy dz 

( 8 ) 

a = A(n • w) 


H = A<|> 


w = (p 2 u,v, w) = (p 2 <t>*,<t>y,<t> z ) = V<J> 


In the above equations a is the source strength and p 
is the doublet strength at any point Q, on the surface S, 
which in our case will be all (for subsonic flow) or part 
(for supersonic flow) of the aircraft surface and wakes. 
These strengths are equal to jumps (discontinuities) across 
the panels of certain flow properties. The source strength 
equals the jump in the normal component of the mass-flux 
vector w. The doublet strength equals the jump in 
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potential, and the gradient of the doublet strength equals 
the jump in tangential component of velocity. The values 
of these strengths are the (as yet unknown) amplitudes of 
the source and doublet singularity solutions appearing in 
equations (3) and (4). Here, these source and doublet 
solutions, when multiplied by a constant k _1 , are denoted 
as K<j and K^, respectively (K is used to denote that the 
singularities are called kernels). For Mo» < 1, k = 4 jc, and 
S is the entire surface of the aircraft and wake(s). For 
M„o > 1, k = 2 jc, and S is that portion of the aircraft 
surface and wake(s) that lies in the upstream Mach cone 
emanating from the influenced point P. 

Equation (7) shows us that the velocity potential at 
point P is related to the source and doublet distributions 
on S, and to the spatial distribution of V 2 <]> in the volume 
V bounded (wetted by) both sides of S. If <t>p is con- 
structed according to the surface integral terms in equa- 
tion (7), that is, 

<|,p = JJjaK a + pKjdS Q (9) 

then, because equation (7) is an identity, it follows that 

JJJ v (v 2 q>) K a d v -° (10) 

Since K<j is a function of the arbitrary point P, V <)> must 
be zero. Thus, construction of 4>p according to equa- 
tion (9) implies that equations (1) and (2), the Prandtl- 
Glauert equation, has been satisfied throughout V. 


°(*j) = <y j(^ T )) = X i 

(12a) 

F(*j)=l l j(^ T l) = x ? 

(12b) 


where the unknown constants Xj and Xj* are called source 
and doublet singularity parameters, respectively, for panel 
j, and (£,q) are local surface coordinates associated with 
the panel. 

Once the functional form for Oj(£/n) and ltj(^t|) are 
specified, equations (9) and (11) can be integrated over 
each panel (a nontrivial task) so that <(>p and v P are 
expressions involving only the unknown singularity 
parameters. If P is made a control point (a panel point at 
which a boundary condition will be imposed) with index i, 
equations (9) and (11) give the potential and velocity at 
that point in terms of (as influenced by) the source and 
doublet distributions of the single panel j (see fig. 7). Note 
that the fixed point P and the variable point Q of the ana- 
lytic formulation correspond to control point i and panel j, 
respectively, in the discretized implementation. Summing 
the effects from all the panels on the aircraft surface gives 
the potential and velocity at control point i in terms of the 
total number (N) of singularity parameters. 

If Vy denotes the velocity at control point i, owing to 
the source-doublet distributions at panel j, then the veloc- 
ity at point i owing to all N panels is 

N 

03) 

j=l 


Equation (9) is the basic starting formula for panel 
methods using sources and doublets. If the source and 
doublet strength distribution is known (we will see how 
this is done later), then the velocity at point P is obtained 
from equation (9) by differentiating with respect to the 
coordinates at P, that is 

v P = V P (Op = JJ[aVpK a + uVpK,]dSQ (11) 


If the panel associated with control point i is to be a solid 
(impermeable) panel represented by a zero normal com- 
ponent of the total velocity, then the boundary condition is 


V r ni = 


f N ^ 

i =1 


fij = 0 


(14) 


Equations (9) and (1 1) are used to generate influence- 
coefficient equations that relate source and doublet 
strengths at particular points Q on the surface S to the 
potential and velocity at field points P. The basic idea is to 
break S into a collection of panels I and to assume a 
functional form for a and p over each panel. For example, 
a constant-strength source-doublet panel with index j is 
given simply by 


Thus, for control point i, we have 
N 

= 05) 

j=l 

and since the v^ are known (from eq. (1 1)) in terms of the 
N singularity parameters, equation (15) can ultimately be 
expressed as the single equation 
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N 

^VIC ij Xj = b i (16) 

j=l 

where the VICy are called velocity influence coefficients 
and b^-V^fq 

Repeating the process for a total of N control points 
and applying boundary conditions at each of these points 
leads to the equation 

[AIC]{X} = {b} (17) 

In this equation, [AIC] is called the matrix of aerody- 
namic influence-coefficients, {X} is the vector of 
unknown singularity parameters, and the elements of {b}, 
the so-called right-hand-side quantities, are known from 
the boundary conditions. Solving for the singularity 
parameters then makes it possible to compute (from 
eqs. (9) and (11)) the potential and velocity distributions, 
and hence the pressures acting at the paneled approxi- 
mation to the surface S. Integrating the pressures and their 
moments yields the resultant aerodynamic force and 
moment. 

To summarize the basic ideas, a panel method uses 
the fact that sources and doublets are solutions to the 
Prandtl-Glauert equation. The numerical procedure is as 
follows: 

1. Break the aircraft surface into an assemblage of 
panels 

Z Use the panels to create source-doublet distribu- 
tions in terms of singularity parameters {X}, whose values 
are to be determined 

3. Form influence-coefficient expressions for the 
potential or velocity or both at each of N panel control 
points owing to the source-doublet distribution of any 
panel; for supersonic flow the zones of silence must be 
accounted for 

4. Use the influence-coefficient expressions to 
enforce boundary conditions at N points, giving the N x N 
system of equations 

[AIC]{X} = {b} 

5. Solve for the singularity parameters {X}, then 
compute the potential and velocity at any point of interest 


6. Knowing the velocity, compute the panel pres- 
sure distributions 

7. Integrate the pressure distributions to obtain 
forces and moments 

The above process, in one form or another, is common to 
all panel methods. Selecting a way to break the surface S 
into panels X, choosing and implementing functional 
forms for a and |i, and evaluating the influence- 
coefficient surface integrals are the fundamental tasks. 

In section 4 we discuss how various codes differ in 
the selection of X, and of a(£,ii) and p(£,rj). Section 5 
describes various means of specifying boundary condi- 
tions and other aspects of modeling physical problems. 
Appendix A gives additional details about how the basic 
approach described above is actually implemented into 
the PAN AIR code. This will lead to specific surface inte- 
grals that must be evaluated to obtain the influence- 
coefficients. The analytic evaluation of some of these 
integrals for subsonic flow is given in appendix B. 

The reader is warned that the above description has 
ignored some subtle points that will be addressed later. 
Specifically, a panel with both source and doublet distri- 
butions requires two boundary conditions per panel. These 
boundary conditions generally require some statement 
about the flow on each side of the panel. (There is an 
exception in supersonic flow.) The boundary conditions 
must also produce a well-posed mathematical problem, 
that is, one with a unique solution; otherwise the AIC 
matrix will be singular. Also, some of the analytic expres- 
sions for the potential and velocity can become singular 
(blow up) when the influenced point is a point of the 
influencing panel (when i = j); hence use of the name 
source-doublet "singularities.” Boundary conditions will 
be discussed in section 5, and improper integrals (those 
having a nonintegrable singularity) will be discussed in 
appendix C. 

4. DISTINCTIONS BETWEEN PANEL CODES 

Having outlined the basic approach used by panel 
methods, let us now look at various implementations. Dif- 
ferent ways of approximating the actual aircraft surface 
with surface panels will be described first, followed by the 
selection of singularity distributions and their functional 
forms. Finally, a summary of several specific panel codes 
will be presented. 
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Surface Geometry 

Current actual-surface panel codes generally allow 
the user to break the configuration into logical pieces such 
as forebody, canopy, and the wing upper and lower sur- 
faces. In the PAN AIR code, such logical pieces are called 
networks of panels (fig. 8). A component such as a wing 
can also be broken into several networks. For example, 
two different models for the F-16XL wing upper surface 
are shown in figure 9. The 3-network model is appropriate 
if control-surface deflections are not of interest; the 
10-network model allows leading and trailing-edge 
control surface deflections to be made. 

Panel codes typically use arrays of M x N surface 
grid points to define the comer points of panels, as illus- 
trated in figure 10. Since the four panel comer points are 
not, in general, coplanar, a single flat panel cannot be 
used to connect the four comer points. 

The lower-order (constant-strength) codes usually use 
some kind of flat ‘‘average” panel defined by the four 
comer points. Using flat quadrilateral panels to represent 
a curved surface causes gaps to exist between the edges of 
adjacent panels (fig. 11, courtesy of John Hess). This 
approach begs the natural question about “leaks” through 
the gaps. Actually, for constant-strength panels, the panels 
themselves leak rather badly everywhere except at the 
control points x, where the discrete boundary conditions 
are imposed. The constant-strength, singularity influence- 
coefficient-computed velocities at any points on the panel 
other than the control point are not generally reliable. Sur- 
face and near-surface values of velocity away from the 
control points are usually obtained by interpolation. 

Parabolically curved panels have been used in at least 
three higher-order, subsonic-flow codes (refs. 31-33). For 
a given number of panels on a curved surface, such panels 
generally provide better accuracy than do flat panels. The 
parabolically curved panels are similar to the flat panels in 
that they still leave gaps between adjacent panels. 

As indicated in the introduction, the difficulties 
encountered with extending panel methods to supersonic 
flow were eventually traced to using doublet-strength dis- 
tributions that were not continuous at panel edges. The 
task of building continuous doublet-strength distributions 
requires that the panel edges themselves be contiguous 
(have no gaps) since the doublet strengths are referenced 
to the panel geometry. 

There are several ways to connect an M x N array of 
arbitrary points witlrfM - 1) x (N - 1) quadrilateral pan- 
els that have no gaps along adjacent panel edges. Probably 


the simplest is to use several piecewise flat subpanels to 
connect four comer points. The most recent version of the 
QUADPAN code uses four triangular subpanels 
(fig. 12(a)), whereas the PAN AIR, HISSS, AND 
MCAERO codes use four triangular subpanels and a flat 
interior parallelogram (fig. 12(b)). Although not shown in 
figure 12(b), the flat parallelogram portion of the panel is 
further subdivided into four more triangular panels 
defined by the diagonals of the parallelogram. This is 
done to facilitate the construction of the continuous 
quadratic doublet-singularity distributions within the 
panel. 

Two codes, SOUSSA and the most recent version of 
the LEV code (leading-edge vortex code, also called the 
free vortex sheet code, ref. 27), use hyperbolic-paraboloid 
panels, (fig. 12(c)). The curved shape of these panels is 
the same as for a structural plate that is loaded as in fig- 
ure 12(d), that is, a twisted shape is produced (ref. 34). 
The side edges remain straight and hence produce no 
gaps. 

Use of the contiguous piecewise flat panels has 
enabled the influence-coefficient-equation surface inte- 
grals (eqs. (9) and (11)) to be evaluated analytically for 
both subsonic and supersonic free-stream flow. Use of the 
hyperbolic-paraboloid panels has been restricted to sub- 
sonic flow (no one has been able to evaluate the surface 
integrals when using these panels for supersonic flow). 

Singularity Distributions 

Most codes use a combination of source and doublet 
distributions on the panels. The primary exceptions are 
the Woodward-Carmichael and USSAERO (also called 
Woodward-2) codes, which use elementary horseshoe 
vortices instead of doublets. The elementary horseshoe 
vortex potential is obtained by integrating the doublet 
potential in the x (streamwise) direction (ref. 35, p. 87). 
The strength of the vortex distribution, also called the sur- 
face vorticity vector y, is related to the doublet strength as 
follows (ref. 2, eq. (A- 3); ref. 1 1, eq. (B.3.9)). 

7(S.n) = n® V|i(S,ti) (18) 

Both y and Vp are in the plane of the panel, n is a unit 
vector normal to the panel, and ® is the cross-product 
operator. 

Source, doublet, and vorticity surface distributions 
cause certain flow properties to be discontinuous with 
respect to the panel surface; that is, there is a jump in the 
flow properties. In general, this means that the velocity 
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vectors on opposite sides of a panel are different. This is 
illustrated in figure 13, where the subscripts U and L refer 
to the “upper” and “lower” sides of the panel, respec- 
tively. For incompressible flow, a source-only panel 
causes the normal component of the velocity at point 
of the panel to jump by an amount given by 

(Vu-V L )-n = o(U) 09) 

and a doublet-only panel causes both the potential and the 
tangential component of velocity to jump. These jumps 
are given by 

<t>u- < t>L = H(^ T l) (2°) 

Vtu-V Tl = V^,ti) (21) 

The jump property given by equation (21) is derived in 
appendix B. Each of the jump properties is derived in ref- 
erence 33. 

The velocity jump across a panel for incompressible 
flow is then given by 

Vu - V L = <rn + Vp. (22) 

The generalization of equation (22) for compressible flow 
is given by equation B.3,29 of reference 12 (version 3.0, 
1990). 

Equations (18) and (22) provide some guidance for 
selecting consistent functional forms for the singularity 
distributions. A constant-strength source panel combined 
with a constant-strength vorticity panel or a linear- 
strength doublet panel will produce a velocity-jump dis- 
tribution that is constant over the panel surface. A linear 
velocity-jump distribution requires a linear-source and 
linear-vorticity or quadratic -doublet distribution. As will 
be discussed later, popular sets of boundary conditions are 
those that include setting the perturbation potential to zero 
on the interior side of panels (the side not wetted by the 
external flow field). Consequently there is no interior 
perturbation velocity, and the velocity jump across the 
panel equals the perturbation velocity on the exterior side. 

As discussed in reference 32, consistency also 
depends on the surface shape used for the panels. Flat 
panels are consistent with constant sources, and parabolic 
panels are consistent with linear sources. In this sense, the 
methods of references 32 and 33 (which use parabolic 
panels, linear sources, and linear vorticity or quadratic 


doublets) are among the few totally consistent higher- 
order formulations. 

Continuity of Doublet Strength 

The reason discontinuous doublet strength can cause 
disastrous numerical problems in supersonic flow is illus- 
trated in figure 14. For a problem simple enough to solve 
in closed form, the actual doublet strength might look 
something like that in the top part of the figure, that is, |i 
varies in a continuous manner. A panel method that does 
not enforce doublet-strength continuity at panel edges will 
produce a solution like that shown in the lower part of the 
figure, that is, the doublet strength jumps at the panel 
edge. Now consider either of the panels separately. In this 
view the doublet strength has a nonzero value at the panel 
edge. The velocity field produced by the doublet distribu- 
tion on either panel is given by the second term of equa- 
tion (11). This term can be integrated by parts (appen- 
dix B of ref. 1 1) to obtain 

v p = ® ( VpR" 1 )^ dri 

+ -^Jn(/)Vp(R- 1 )®d/' ( 23 ) 

where the line integral is around the panel perimeter. 

The surface integral in equation (23) involves the sur- 
face vorticity vector y , and is called the regular part of the 
doublet velocity. The line integral involves the doublet 
strength at the panel edges and is called the line-vortex 
term. The two line- vortex terms from the common edges 
of the two panels in figure 14 produce a single line- vortex 
of strength T(i) = A p,(f), where A p,(/) is the jump in the 
doublet strength along the panel edge. The velocity field 
from this line vortex is spurious since the actual doublet 
strength does not jump at the panel edge; that is, the dis- 
cretization has introduced velocities that should not exist 

In subsonic flow these spurious velocities decay 
rapidly with distance from the edge and usually do not 
cause serious problems. In supersonic flow, these veloci- 
ties persist, their effect propagaung down the Mach cones. 
Consequently, erroneous incremental flows continue to 
exist at control points within the domain of influence of 
the disturbance points, thereby introducing errors in the 
AIC matrix. These errors are frequently serious enough to 
produce a totally incorrect solution for the flow. 

Now, if it were known in advance that the discretized 
doublet strength would always be continuous in value. 
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then the resultant line-vortex term along common panel 
edges of adjacent panels would always be zero. Hence, 
the line-vortex term in equation (23) would not even have 
to be computed, but could be simply discarded. This 
approach is taken in the subsonic-supersonic-flow PAN 
AIR and HISSS codes: the doublet strengths are made 
continuous by construction and the line-vortex terms are 
thrown away. 

A comparison of results for a 3%-thick swept wing at 
Moo = 2.05, based on discontinuous and continuous 
quadratic -doublet distributions, is shown in figure 15 
(figs. 9 and 10 in ref. 36; figs. 48 and 62 in ref. 37). In the 
discontinuous-distribution case, there is a pressure spike 
at about 80% chord of the station just inboard of the tip. 
This is caused by a doublet-strength discontinuity near the 
leading-edge region of the wing tip, the effect of which 
propagates down the Mach cone originating at the tip 
leading edge. The spike is not present in the continuous 
doublet-distribution case. 

Summary of Methods 

Table 1 summarizes the basic features of several 
panel codes. Codes that handle only subsonic flow are 
listed separately from those that treat both the subsonic 
and supersonic cases. The QUADPAN code is listed twice 
because it was originally a subsonic -only code and was 
later revised to do supersonic flow. The year given next to 
each code’s name is the approximate year the code was 
introduced. The panel geometry and the singularity types 
and their spatial distributions are listed for each code, A 
footnote indicates that geometry or singularity type is con- 
tinuous from panel to panel. The table includes what are 
probably the best known, or most generally available 
codes; other codes (often proprietary) also exist. Refer- 
ences 39 through 41, as well as other sources (material 
distributed at the 1985 AIAA Workshop on Aerodynamic 
Analysis Using Panel Methods), compare results from 
several of the codes for various configurations. 

Codes for subsonic-only flow- The lower-order 
Hess code (also known as the Douglas-Neumann code) is 
considered by many to be the first practical implementa- 
tion of the panel method for quite general geometry 
(ref. 2). It actually exists in several versions, the first 
being a source-only version that did not treat lifting 
problems. The lifting case was later added. The Hess 
family of codes is one of the few that enforce the Kutta 
condition with a nonlinear pressure rule. As shown in 
reference 40 this can sometimes give a better trailing -edge 
pressure result than is possible with a linear approxi- 
mation to the pressure rule. 


MCAERO uses higher-order continuous doublet dis- 
tributions on piecewise flat continuous panels (ref. 17). 
One of its important capabilities is the use of analytically 
differentiated (with respect to panel coordinates) 
influence-coefficients. These are used to efficiently 
implement the so-called design problem, that is, the prob- 
lem of determining the geometric shape required to pro- 
duce a specified pressure distribution. 

SOUSS A stands for steady, oscillatory, and unsteady 
subsonic and supersonic aerodynamics (ref. 8). Although 
it produced some supersonic results by using a small but 
nonzero reduced frequency, it turned out to be incapable 
of doing steady supersonic flow. One of its legacies was 
to popularize the so-called Morino boundary-condition 
formulation, which is discussed in section 5. This formu- 
lation led to renewed interest in lower-order subsonic 
codes and resulted in the development of the VSAERO 
and QUADPAN codes. 

VSAERO (from vortex separation aerodynamics 
analysis) is one of the few codes that contain a procedure 
for calculating the shape and location of the trailing-wake 
system (ref. 18). 

The LEV code (leading-edge vortex; ref. 27) is also 
known as the free-vortex sheet code. It was designed 
specifically to model the vorticity shed from sharp leading 
edges of swept wings (fig. 5). 

The higher-order Hess code (ref. 32) is the most 
recent version of the Hess family of codes. It uses 
parabolically curved panels in conjunction with linear 
source and linear vorticity (quadratic doublet) 
distributions. 

Codes for subsonic or supersonic flow- The 
Woodward-Carmichael code (refs. 4, 5) (also known as 
the constant-pressure panel code and the Woodward-I 
code) was described in the introduction. It is still used for 
simple configurations that can be approximated with the 
mean-surface representation of figure 1(b). 

USSAERO (unified subsonic/supersonic aerodynam- 
ics) is also known as the Woodward-II code (ref. 7). It 
differs from the Woodward-I code in that the line distri- 
butions used for bodies were changed to constant-strength 
source panels, and the constant-strength vorticity panels 
were changed to distributions wherein the strength varies 
linearly in the chordwise direction and is constant in the 
spanwise direction. The corresponding doublet strength 
(quadratic-linear) is not continuous, so flat mean-surface 
models of lifting surfaces are usually required for super- 
sonic flow. 
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PAN AIR (panel aerodynamics) is generally consid- 
ered to be the first actual-surface paneling code with reli- 
able numerics for supersonic flow (refs. 9-14). In addition 
to using continuous doublet distributions, it was also nec- 
essary to incorporate forward-weighted splines for the 
supersonic case. This is somewhat analogous to the use of 
forward-differencing in the finite-difference computa- 
tional fluid dynamics (CFD) codes. The basic technology 
in PAN AIR is also used in the HISSS code (higher-order 
subsonic/supersonic singularity method, refs. 16 and 42). 

QUADPAN (quadrilateral panel aerodynamics pro- 
gram) started out as a subsonic-only code and used 
constant-strength sources and doublets (ref. 19). Later, the 
doublets were changed to a linear-continuous distribution 
so that supersonic flow could be handled (from 1985 
AIAA Workshop on Aerodynamic Analysis Using Panel 
Methods). 

It has been well established that continuous doublet 
distributions are essential for reliable numerics in super- 
sonic flow. However, it is probably fair to say that there is 
not a unanimous agreement between panel code develop- 
ers about the need for higher-order approaches for the 
subsonic-flow codes. The advantages that lower-order 
codes offer over higher-order codes are (1) less work to 
derive the influence-coefficient equations; (2) simpler 
coding implementation (because the higher-order 
approaches must relate information involving different 
panels, which leads to special cases and logic); and (3) far 
fewer arithmetic operations, hence lower run costs. One 
reason some subsonic-only codes use higher-order 
approaches is that their developers believe the numerics 
are more reliable than those of the lower-order approaches 
for highly complex geometry. In the case of two codes 
that use parabolic panels with linear sources and quadratic 
doublets, numerical calculations have demonstrated that 
fewer panels are required than for lower-order methods to 
obtain a given accuracy (refs. 32, 33). This seems to be 
especially true for internal flows. 

5. MODELING 

Modeling refers to particular techniques used to 
simulate flow about an object The modeling tools at the 
disposal of a panel-method user are (1) the geometric 
generality that panel codes provide, (2) the use of sources 
and doublets, individually or in combination, and 
(3) boundary conditions. In this section we will discuss 
boundary conditions, their interplay with the source- 
doublet jump properties introduced in section 4, and how 
the combination can be used to model in different ways. 
There is also a brief look at wake modeling. 


In section 2 we discussed the general limitations 
inherent in the Prandtl-Glauert equation derivation, 
namely, the flow is represented as being inviscid, irrota- 
tional (potential), and linear. If the free-stream Mach 
number is supersonic, there are additional geometric 
restrictions; these are illustrated in figure 16. 

In supersonic flow, the Prandtl-Glauert equation 
admits solutions for solid surfaces only if the surfaces are 
swept back more sharply than the Mach cone. Thus, the 
higher the Mach number, the more streamlined must be 
the aircraft. This restriction means that forebodies must be 
pointed, not blunt Wing leading edges can be blunt only 
if they are swept behind the Mach cone (a so-called sub- 
sonic leading edge). 

Surfaces can be swept at angles smaller than the 
Mach cone only if they do not represent solid surfaces. 
These so-called superinclined panels (ref. 1 1) are used for 
nacelle inlet faces and nozzle exit planes. Superinclined 
panels use both sources and doublets, and both boundary 
conditions must be prescribed on the downstream side of 
the panel. Numerical experiments indicate that these 
downstream boundary conditions must also be of a spe- 
cific type, namely, that the potential and its normal deriva- 
tive must be specified. This is analogous to an initial 
value problem. Superinclined panels have no upstream 
influence. The panels generate a downstream flow only 
and simply absorb any flows that run into them. Refer- 
ence 11, appendixes A and B in reference 12, and refer- 
ence 37 are recommended sources of information on the 
supersonic aspects of modeling. 

Mass-Flux and Velocity Boundary Conditions 

The physical description of a real flow at a surface is 
given by the no-slip boundary condition 

V = 0 (24a) 

or 

(pv) = 0 (24b) 

These equations state that the total velocity vector V, or 
the mass flux vector (pV), is zero at a solid surface. For 
inviscid flow, the tangential component of the velocity 
cannot be prescribed (unless the pressure is known) and 
equations (24) are replaced with the zero normal-flow 
boundary conditions, 

V • n = 0 (25a) 
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or 


V'n' = 0 


(30) 


(pv)n = 0 (25b) 

which must be supplemented with a Kutta condition for 
subsonic trailing edges. 

Equation (25b) is nonlinear since the density is a 
function of the unknown velocity; a linear approximation 
to equation (25b) is (from sec. 1. 1 1 of ref. 30) 

pV^pJ^ + p^V* (26) 

Equation (26) involves neglecting terms that are of the 
same order of magnitude as those neglected in deriving 
the Prandtl-Glauert equation. Dividing by poo gives 

W = V = Vo, + w (27a) 

P~ 

where, from equation (8) 

W = V<J> = (p 2 <l> x .<> y ,«t» z ) (27b) 

The quantities W and w are called the linearized total 
and perturbation mass flux, respectively. Thus, equation 
(25b) is approximated by the mass-flux boundary 
condition 

W*n = 0 (28) 

which is linear in the components of the perturbation 
velocity, and depends on the free-stream Mach number. If 
the panel is not meant to be impermeable, for example, at 
an inlet face, then the right-hand side of equation (28) is 
replaced with a specified nonzero value. 

For a fixed subsonic Mach number, a solution for 
<|>(x,y,z) can be obtained in either of two ways. The first is 
to solve the “true” compressible problem, that is, solve the 
Prandtl-Glauert equation as it stands, using the mass-flux 
boundary condition applied to the true geometry. The sec- 
ond (and better known) way is to solve the equivalent 
incompressible problem by using the Prandtl-Glauert 
transformation to convert equation (1) to Laplace’s equa- 
tion (sec. 7-1 in ref. 35): 

4>X V + 4>y'y' + <t>zV = 0 (29) 

using the velocity boundary condition 


applied to the transformed geometry. In equations (29) 
and (30) the primes indicate that all the variables 
are in the transformed coordinates, for example, 
V' = + (<K',<l>y'»<l>z') • Transforming the solution for <j>' 

back to the true (physical) variables gives the solution for 
4> for the true geometry. The important point is that the 
solution for 4> obtained from the first approach will be the 
same as that obtained from the second approach (ref. 43). 
Thus, solving for the flow about the true geometry by 
using the subsonic Prandtl-Glauert equation with mass 
flux (instead of velocity) boundary conditions, is mathe- 
matically the same as solving the equivalent incompress- 
ible problem with velocity boundary conditions applied to 
the transformed geometry. The first approach has the 
advantage that it can also be used for supersonic flow 
(where there is no equivalent incompressible problem). 

Note that the linearized mass flux defined by equa- 
tions (27) is not parallel to the true mass flux pV, and, 
hence, is not parallel to the velocity vector V = V*, + v. 
This inconsistency is one of the prices paid for the lin- 
earization. The consequences of this will be illustrated 
with the following two examples. (Theoretical discussions 
of mass-flux and velocity boundary conditions appear in 
refs. 44A6.) 

The first example demonstrates the accuracy of the 
mass-flow boundary condition and the Prandtl-Glauert 
equation for supersonic flow over a thin wedge. It also 
demonstrates the jump properties of the mass-flux and 
velocity vectors across shocks predicted with the Prandtl- 
Glauert equation. 

Panel-method results obtained from the PAN AIR 
code are shown in figure 17 (from ref. 9). The boundary 
condition applied to the exterior side of the panels is the 
linearized mass flux condition W • n = 0. This forces the 
linearized mass flux vector W , instead of the resultant 
velocity vector V, to be parallel to the wedge faces. Pres- 
sures have been computed with the isentropic and the 
second-order pressure rules (eqs. (8.10) and (8.11), 
respectively, in ref. 23). Also shown are the pressure 
distributions computed from shock -expansion theory and 
from classic linear thin-airfoil theory (secs. 4.16 and 4.17, 
respectively, in ref. 23). The PAN AIR results are very 
close to those of nonlinear shock-expansion theory, 
having a greater pressure magnitude on the frontward- 
facing compression surface than on the rearward- facing 
expansion surface. This contrasts with the classic linear 
thin-airfoil theory, which predicts equal and opposite 
pressures on the two inclined faces. 
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How does one explain that PAN AIR, which solves 
the linear Prandtl-Glauert equation, produces results that 
agree more closely with nonlinear shock-expansion theory 
than with classic linear th in-airfoil theory? Apparently, for 
this particular problem, the velocities are small enough to 
make negligible the nonlinear terms in the differential 
equation. This conclusion is justified by the close agree- 
ment between the isentropic and second-order pressure 
rules. (It is good practice to always compute the pressures 
by these two rules; a substantial disagreement in the 
results is a sure sign that the small-perturbation assump- 
tions of the Prandtl-Glauert equation are being violated.) 
Since classic thin-airfoil theory is also based on the 
Prandtl-Glauert equation, but in addition depends on 
complete linearization of the boundary conditions, the 
explanation must reside in the boundary conditions. 
Although PAN AIR uses^the linearized mass-flow W in 
the boundary condition W • n = 0, the unit normals n act 
at the actual wedge surfaces. In the classic linear-theory 
solution the boundary conditions are not applied to the 
true geometry. Instead, the boundary conditions are 
applied along the straight upper side of the wedge, and the 
flow is made to turn through the wedge angles, as in the 
mean-surface modeling of the Woodward-Carmichael 
code (see sec. 1). So, this appears to be an example of an 
instance in which the solution is more affected by approx- 
imations made to the boundary conditions than it is by 
those made to the differential equation. 

The mass-flux and velocity vectors at two points on 
opposite sides of the Mach line emanating from the wedge 
apex are shown in figure 18. We first note that for expan- 
sive flow about a comer, such as at the wedge apex, the 
Prandtl-Mayer solution predicts an expansion fan. That is, 
away from the apex the flow properties change gradually 
through the fan. However, the Prandtl-Glauert equation 
predicts abrupt changes that occur along the apex Mach 
line. This ‘‘expansion shock” approximation to an expan- 
sion fan is a poor representation at large distances from 
the apex, but is accurate near the apex. Our point here is 
to illustrate what happens to a panel-method solution 
across such a shock (whether expansive or compressive). 
Both the mass-flux and velocity vectors are discontinuous 
across the Mach line; they jump in both magnitude and 
direction. For the mass flux, the tangential component 
jumps, and the normal component is continuous. For the 
velocity, the tangential component is continuous, and the 
normal component jumps. This is the same behavior as 
predicted by the nonlinear Rankine-Hugoniot relations, 
thus providing some measure of confidence in the lin- 
earized mass-flux approach. 

The second example (fig. 19) shows the pressure at 
any point on 10° and 15° half-angle cones at zero angle of 


attack, for free-stream Mach numbers between 1.0 and 4.0 
(fig. 4.1 in ref. 15). At the lower Mach numbers, the mass- 
flux-boundary-condition solution agrees more closely 
with the exact Euler-equation solution than does the 
velocity-boundary-condition solution. As the Mach num- 
ber increases, the mass-flux and velocity-boundary- 
condition solutions become less accurate. The mass-flux- 
boundary -condition solution rapidly diverges from the 
exact solution, and crosses the less rapidly diverging 
velocity-boundary-condition solution. As the cone angle 
increases, the Mach number at which the mass-flux and 
velocity-boundary-condition solutions cross one another 
becomes smaller. This supersonic “thick-body” behavior 
is responsible for the fact that when panel methods are 
applied to nonslender fighter forebodies, the velocity- 
boundary condition often gives superior answers to the 
mass-flux boundary condition in the region just ahead of 
the canopy (refs. 47, 48). For such cases, the mass-flux 
boundary condition can actually produce negative pres- 
sure coefficients, as suggested by the 15° cone solution 
behavior of figure 19. For bodies that are adequately slen- 
der, the linear Prandtl-Glauert equation with mass-flux 
boundary conditions can provide good answers, for exam- 
ple, the Mach 1.6 pressures on the B-l forebody presented 
in references 36 and 37. The important conclusion to be 
drawn from figure 19, and from other examples in ref. 15, 
is that inviscid supersonic flow solutions based on the 
Prandtl-Glauert equation can be substantially in error if 
the Mach number, thickness, or angle of attack are too 
large. In such cases, codes based on nonlinear theories 
(for example, the full -potential equation or Euler equa- 
tions) must be used for reliable answers. 

Interior Potential 

An actual-surface panel model of an aircraft generally 
produces a set of panels that separates space into two or 
more distinct regions: enclosed interior volumes and an 
external volume extending from infinity to the external 
side of the panels. The flow in the external volume corre- 
sponds to the physical flow field being modeled. The flow 
in the interior volumes is fictitious but, as will be shown, 
can be used to advantage. 

Newcomers to panel methods often find the idea of 
an internal flow field to be strange in that no such flow 
exists inside a real wing. However, it must be remembered 
that we are using sources, doublets, and boundary 
conditions to create the flow fields, and flow will exist on 
either side of the source-doublet panels. The external- 
internal flow fields are in general independent of one 
another; they depend on the boundary conditions on the 
external-internal sides of the panels. These include direct 
boundary conditions explicitly imposed on each side of 
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the panel and indirect boundary conditions resulting from 
some specification of the source or doublet strength. The 
indirect boundary conditions arise because of the jump in 
flow properties produced by sources and doublets. 

The boundary condition implied by the equation (20) 
jump property is the cause of a common mistake made by 
new users of panel methods, namely, not specifying the 
potential at a point of an interior volume. When an interior 
volume is totally enclosed by panels that all have doublet 
distributions, the potential must be specified at one or 
more interior points, and this must be done with a source 
distribution. If this is not done, the problem is ill-posed, 
that is, it does not have a unique solution, and the AIC 
matrix will theoretically be singular (will not possess an 
inverse). (In practice, numerical round-off error often 
changes a theoretically singular matrix to just an ill- 
conditioned matrix.) 

The reason for the problem being ill-posed is gener- 
ally explained as follows. At x = °° the value of the poten- 
tial owing to a source or doublet is zero (this is called the 
far-field boundary condition). At the exterior side of pan- 
els enclosing an interior volume the potential will have 
some distribution <|>ext( x *y> z ) that depends on the surface 
boundary conditions and on the resulting source-doublet 
strength distribution. On the interior side of the panels, the 
potential is ^intO^y* 2 )* which by equation (20), differs 
from the external distribution by an amount equal to the 
doublet-strength distribution p(J;,T|). Now, if at some inte- 
rior point, a constant c, arbitrary except that it must satisfy 
any boundary condition associated with the point, were to 
be added to (fentO^y* 2 )* then fanfayz) + c would also sat- 
isfy the Prandtl-Glauert equation and, hence, be another 
solution. The way to make the interior solution unique is 
to specify a value for the potential at an interior point, 
thus determining the constant c. (The same argument can 
be made for the exterior side of the panels, which requires 
that the arbitrary constant be zero to satisfy the far-field 
boundary condition tyo* = 0. Thus, no explicit specification 
of the potential at an exterior point has to be made.) 

The interior potential at a point must be specified 
with a source panel, and not a doublet panel, because of 
the “physically” different behavior between sources (or 
sinks) and doublets. With any paneled geometry, the 
boundary conditions are enforced only at the control 
points. Hence, when a finite number of control points is 
used to impose boundary conditions there will be some 
degree of "leakage.” That is, the integral of the source 
strength (jump in normal component of mass flux) over 
the surface enclosing the interior volume will not be 
exactly zero. When a source panel is used to specify 
velocity potential in the interior domain, the panel is 


capable of generating or consuming fluid to conserve 
mass that is leaked out of or into the interior domain. (The 
source strength is one of the unknowns that are solved 
for.) The doublet panel does not have the ability to gen- 
erate or swallow fluid (at least, in a net sense), and conse- 
quently is not capable of handling the leaked flow. 
Numerical examples demonstrating this behavior are 
given in cases 5 and 6 of reference 49 (pp. 6-7). 

Another interesting example in reference 49 is case 7. 
Here, a thick wing is modeled using doublets-only for the 
wing upper and lower surfaces and sources-only for the 
panels closing off the wing-tip opening. Using only a 
single type of distribution on all panels allows only one 
boundary condition per panel, which was chosen to be 
W ■ n = 0 on the exterior panel sides. The interior poten- 
tial was not specified anywhere inside the wing. At first 
glance this also appears to be an ill-posed problem, but it 
is not The source-only tip-closure panels do not produce 
a potential jump. Therefore, as far as the potential is con- 
cerned, this model behaves as if there were no tip closure, 
and in this sense the wing interior is not a separate domain 
from the exterior domain. Hence, the constant <j><x> = 0 
applies to both the exterior and interior part of the single 
domain. The zero-normal-flow boundary condition 
applied at the tip, however, does separate the flow into 
two regions. 

In actual practice, the interior potential is set not at 
just a single point, but at the interior-side control point 
locations of all panels. If the boundary condition 
<t>int( x »y» z ) = 0 is used at all such control points there will 
be no internal perturbation flow (assuming an adequate 
number of panels), giving uniform free-stream flow in the 
interior volume. If both sources and doublets are used 
with each panel, then another boundary condition is 
needed, and it is selected to control the flow on the exte- 
rior sides of the panels. Two alternative ways of doing 
this, for flow about a solid (impermeable) aircraft surface, 
are described next The first is a direct approach, the sec- 
ond is an indirect one. 

The direct approach is illustrated in figure 20(a). On 
the exterior panel side, the zero normal mass flux bound- 
ary condition (eq. (28)) is imposed. The resulting bound- 
ary condition pair at each panel control point is 

(Wn) =0 (31a) 

V /ext 

Oint = 0 (31b) 
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The indirect approach (sometimes referred to as the 
Morino formulation) is illustrated in figure 20(b). Here, 
equation (31a) is replaced by a condition imposed on the 
source strength, and the boundary condition pair is 

a = - • n (32a) 

4>int = 0 (32b) 

The equivalence between these two boundary condition 
pairs can be shown by starting with the definition of the 
source strength (see eq. (8)): 

® = (wext-w int )n (33) 

If the perturbation potential is zero everywhere in the 
interior region, then its derivative is zero in every direc- 
tion and is then zero. (Equation (32b), applied to the 
interior side of every panel normally produces a good 
approximation to this condition.) Then equation (33) 
becomes (using the definition for W given by eq. (27a)) 

a = wjj-n 
= (Wu-V«)n 

= W u n-V.n (34) 

Now imposing the W u n = 0 condition sought 
(eq. (31a)), we get equation (32a). 

The indirect approach has two features that reduce 
run cost. First, only the influence-coefficients for the 
potential, and not the three components of velocity, have 
to be computed. Second, since the source strengths are 
specified, only the doublet strengths have to be solved for, 
cutting the size of the [AIC] matrix roughly in half. (The 
reason the size is not necessarily cut exactly in half is that 
some codes use more doublet unknowns than source 
unknowns, e.g,, for wakes, and for doublet matching in 
the continuous-doublet codes.) 

The direct approach requires that both potential and 
velocity influence-coefficients be computed, and both 
source and doublet strength must be solved for. For a 
given number of panels this causes the direct approach to 
be more costly than the indirect approach. The advantage 
of the direct approach is that it is sometimes more accu- 
rate than the indirect approach. This is apparently because 
the indirect approach depends on the internal perturbation 
potential being zero everywhere for equations (32) to 
accurately represent equations (31). In practice, the poten- 


tial is set to zero only at discrete interior points; conse- 
quently, the indirect formulation sometimes produces 
more error than does the direct approach. 

An additional advantage, common to both formula- 
tions, is the property that the internal flow is everywhere 
uniform (it equals the free-stream velocity). For super- 
sonic free-stream Mach numbers this prevents the forma- 
tion of internal disturbances that otherwise would propa- 
gate along Mach lines and reflect from the internal sides 
of the panels. If not eliminated, these reflections can cause 
severe internal flow disturbances that are “felt through the 
panels” and degrade the external flow-field solution. 
Although there are no Mach-line disturbances in subsonic 
flow, large velocities are normally produced by the line- 
vortex behavior of constant-strength doublet panels. 
These large velocities are eliminated in the internal flow 
region by the internal potential being set to zero. This is 
apparently why solutions obtained with the newer 
constant-strength source-doublet formulations are not as 
sensitive to the panel layout as were earlier constant- 
strength codes (which did not use a zero interior- 
perturbation potential). 

A final example uses doublet-only panels to represent 
a wing, including the tip closure. Since the interior vol- 
ume is totally enclosed with doublet panels, the potential 
must be set in the interior volume for the problem to be 
well-posed. Also, having only one type of singularity dis- 
tribution means that only one boundary condition can be 
employed per control point. The direct approach would be 
to set the interior potential with an additional source panel 
located inside the wing, and specify W * n = 0 on the exte- 
rior. A clever indirect approach is to simply set the total 
(instead of the perturbation) potential to zero on the inte- 
rior side of each panel. For sufficiently dense paneling, 
this will make the total potential close to uniformly zero at 
all interior points. Hence the gradient of the total potential 
in any direction, which is the total (linearized) mass-flux 
in that direction, will also be zero. Consequently, the total 
mass-flux component normal to the interior side of the 
panels will be zero. Then, since doublets do not produce a 
jump in normal mass flux, the normal mass flux on the 
exterior sides of the panels will also be zero, that is, 
W • n = 0 on the exterior is produced. 

Wakes 

Wake panels are used to enforce the Kutta condition 
at sharp (usually the trailing) edges of lifting surfaces. 
Since wakes trail downstream from these edges, they also 
influence the flow experienced by downstream compo- 
nents. For example, the load experienced by a wing will 
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depend on the proximity of the wake from an upstream 
canard. 

Figure 21 shows a PAN AIR model of a supersonic 
fighter wing-canard concept (ref. 50). Two of three wake 
models tried for this configuration are shown in figure 22. 
Simple, flat wakes for both the wing and canard are used 
in the first model, the wakes being positioned approxi- 
mately in the wing-canard mean plane (fig. 22(a)). For 
each angle of attack, the resulting span loadings were then 
used by a multiple lifting-surface code to compute the 
approximate rolled-up shape of the canard wake. The 
rolled-up canard wakes were used for the second model 
(fig. 22(b)). The third model (not shown) again used flat 
wakes, but for each angle of attack the canard wake was 
aligned with the free-stream velocity direction. The results 
in reference 50 indicate that for this configuration, the 
second and third models give essentially the same lift and 
moment, and these agree more closely with wind-tunnel 
results than do the results from the first wake model. Lift 
and pitching moment coefficients at Mach 2.2, obtained 
from the third wake model, are shown in figure 23. There 
is a noticeable difference in the PAN AIR-predicted lift- 
curve slopes for the isentropic and second-order pressure 
rules, indicating that nonlinear effects are beginning to 
become important. Also shown are results from the 
USSAERO (version B) program; it does less well at pre- 
dicting the moment. 

Another wing-canard configuration is shown in fig- 
ure 24 (ref. 51). The effect of aligning a flat-canard wake 
model with the free stream instead of with the chord plane 
(unaligned position of fig. 25) for this close -coupled 
wing-canard is illustrated in figure 26. This figure shows 
the spanwise circulation distribution of the canard and the 
wing for the aligned and unaligned flat-canard wakes. The 
canard-wake position does not significantly affect the lift 
distribution along the canard. It does, however, have an 
important effect on the wing-lift distribution. Moving the 
canard wake from the unaligned position to the aligned 
position causes two major changes in the flow over the 
wing, the dominant one of which is an increase in wing 
lift inboard of the canard-tip station. This is due to the 
diminished canard-wake downwash field, raising the 
effective angle of attack of the inboard wing section. The 
secondary effect is a loss in wing lift outboard of the 
canard-tip station. This loss in lift is due to diminished 
spanwise velocity imparted on the upper surface of the 
wing by the canard wake, which is due to the increased 
distance between the canard wake and the wing. The net 
effect of these changes caused by the aligned wake model 
is an increase in total lift over that predicted by the 
unaligned wake model. 


The final example, taken from reference 52, is a PAN 
AIR flaps-down analysis of the Boeing 737-300 (fig. 27). 
The actual and computational flap geometries for the 
flaps- 15 setting, used for most takeoffs, are shown in fig- 
ure 28(a). The assumed flap, wing, and slat wake posi- 
tions are shown in figure 28(b). Leading- and trailing- 
edge spanwise geometry discontinuities owing to the 
deflected flaps and slats also had to be treated; these are 
discussed in reference 52. 

The drag buildup was obtained from a combination of 
methods. The profile drag for everything but the wing was 
estimated from handbook data. The wing profile drag was 
estimated from a two-dimensional multielement-airfoil 
panel code coupled with a viscous model (ref. 53). The 
induced drag was obtained from PAN AIR’s surface- 
pressure integrations. Figure 29(a) shows the drag contri- 
buions for both the flaps- 15 (FI 5) setting, and a lower set- 
ting called flaps-1 (FI). The resulting computational lift- 
to-drag ratios are compared with flight data in fig- 
ure 29(b). 

6. TRANSONIC FLOW 

Panel-method codes have existed for about 25 years 
and are still the only codes routinely used to analyze flow 
about complex three-dimensional configurations. Their 
general inability to solve nonlinear problems is a serious 
drawback, however, and most current CFD research 
involves finite-difference, finite- volume, or finite-element 
approaches to solving nonlinear flow equations. In the 
United States, the finite-difference and finite-volume 
approaches currently seem to be in most favor. Unfortu- 
nately, most finite-difference and finite-volume 
approaches require well-structured flow-field grids that 
conform to the surface of an aircraft. It “is the difficulty in 
generating ‘suitable grids’” (ref. 54) that is the major 
technical obstacle to routinely computing in viscid tran- 
sonic flow about realistic aircraft. In contrast, it is rela- 
tively easy to produce panel-method-type grids that are 
only on the aircraft surface. 

This raises a question: Is it possible to combine the 
surface geometry grid used by panel methods, with some 
easy-to-generate flow-field grid and use the combination 
to solve nonlinear fluid flow problems? The answer is yes, 
as demonstrated by the TRANAIR code and the work that 
led to it (refs. 55-59). This section provides a brief intro- 
duction to the technique, and how it evolved. To the user, 
TRANAIR appears to be a panel code since the input is 
panel-code-like; however, the solution techniques 
TRANAIR uses are not those of a panel code. 
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The basic approach used in TRAN AIR is to embed 
the surface panels in a rectangular box of grid points, as 
shown in figure 30. The initial formulation (sec. II of 
ref. 56) was based on Green’s third theorem (eq. (7)), and 
combined the surface-integral-generated influence-coeffi- 
cients with the volume integrals. The rectangular grid was 
used to evaluate the volume integral, for every point of 
the rectangular and surface grids, with fast Fourier trans- 
forms. This Green’s theorem approach was able to solve 
the nonlinear full-potential equation if there were no 
shocks, but proved to be unstable when shocks were pres- 
ent. By resorting to more powerful mathematical methods 
involving optimization, supercritical results using 
influence-coefficients were subsequently obtained 
(refs. 55, 56); however, the computational cost was 
extremely high. 

The most recent version of TRANAIR (refs. 57-59) 
still uses panels and the rectangular grid, but does not use 
influence-coefficients. Instead, the cells formed by the 
rectangular grid are used to discretize the full-potential 
equation with tri-linear basis-function finite elements. The 
surface panels, which slice through some of the finite- 
element cells, alter the finite-element discretization in the 
vicinity of boundary surfaces. 

At the perimeter of the rectangular grid, the equation 
set being solved changes from the full-potential equation 
to the Prandtl-Glauert equation. As a consequence, the 
rectangular computational grid need only encompass the 
nonlinear flow regions (which are only near the aircraft). 
The far-field boundary condition of zero potential at infin- 
ity is automatically satisfied by the discrete Green’s func- 
tion (for the Prandtl-Glauert equation) used in the formu- 
lation. Consequently, the solution domain extends to 
infinity, even though the computational grid is finite, as 
indicated in figure 30. 

The finite-element discretization yields a set of non- 
linear algebraic equations which are solved iteratively, 
using Newton linearization, multiple preconditioners, and 
an optimization algorithm called GMRES (generalized 
minimal residual). Details are given in references 58-60. 
The computer run cost based on this approach is much 
less than with the original influence-coefficient/ 
optimization approach. 

The input to TRANAIR is essentially that of the PAN 
AIR code, that is, the surface grids of panel comer points 
supplemented by the box of rectangular flow-field grid 
points. This enables transonic flow to be computed about 
very complex configurations without having to generate a 
surface-conforming flow-field grid. For example, 
TRANAIR has been used to compute transonic full- 


potential solutions for the F-16A, using the rectangular 
grid box and paneled geometry of figure 30. For the half 
geometry shown, the grid box contains 129 x 33 x 33 
points, and the aircraft contains about 3,500 panels. Two 
views of the surface paneling are shown in figure 31. In 
this model, the wing-tip missiles and launchers are not 
included. 

F-16A supercritical wing-pressure results (ref. 57) are 
shown in figure 32. The free-stream Mach number is 0.9 
and the angle of attack is 4°. The experimental data for the 
outboard station indicates separated flow near the trailing 
edge. This is probably a result of the wing-tip missiles and 
launcher that were part of the wind-tunnel model. The 
wind-tunnel data indicate a shock at approximately 75% 
chord for the four inboard stations. TRANAIR also indi- 
cates a shock, but it is slightly downstream of the shock 
predicted by the wind tunnel. This result is generally 
expected from a conservative full-potential solution, a 
result of the absence of a boundary-layer correction 
(fig. 40 in ref. 61). The shock predicted by TRANAIR is 
smeared over 5 to 6 grid cells, and at the 2 outboard sta- 
tions where there are only about 11 rectangular-grid 
points spanning the chord, the shock is completely 
washed out. Increasing the grid density in the x-direction 
(streamwise) greatly reduces the shock smearing. To bet- 
ter resolve rapidly varying flow behavior without increas- 
ing the flow-field grid density everywhere, techniques 
such as local grid refinement, wherein individual cells of 
the rectangular grid are subdivided into smaller cells, or 
higher-order finite-element-basis functions are needed. 

The coarseness of the rectangular grid is particularly 
evident in the leading-edge region of the wing. An exam- 
ple is shown in figure 33, at about 70% semispan. The 
first 18% of the wing chord is spanned by only three cells 
of the rectangular grid. Consequently, with this coarse 
uniform grid, TRANAIR fails to capture the leading-edge 
pressure peaks. The subsequent addition of local grid 
refinement has enabled these pressure peaks to be 
resolved (ref. 62). 

7, CONCLUDING REMARKS 

An attempt has been made to give an overview of the 
basics of the panel method and to provide some fairly 
specific details on how the basics can be implemented. 
The tools at the disposal of the panel- method user are 

(1) surface panels of source-doublet- vorticity distributions 
that can represent nearly arbitrary geometry, and 

(2) extremely versatile boundary-condition capabilities 
that can frequently be used for creative modeling. 
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Panel methods have reached a relatively mature stage 
of development. Their fundamental limitation is that they 
solve only linear differential equations. Even so, they are 
widely used in the aerospace industry because they can be 
used to model extremely complicated geometry. 

Recently, panel-method technology has been com- 
bined with other procedures to solve transonic flow prob- 
lems. For example, the TRANAIR code has been able to 
solve the full-potential equation for the F-16A at super- 
critical Mach numbers. It does this by combining surface 


panels with a rectangular flow-field grid, thereby 
eliminating the often difficult task of creating three- 
dimensional, surface-fitted, flow-field grids. To the user, 
TRANAIR appears to be a panel code because the input is 
essentially the same as that of a panel code; however, the 
actual solution process is based on finite -element and 
optimization techniques. 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, CA 94035-1000, February 13, 1990 
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APPENDIX A 


PAN AIR IMPLEMENTATION OF SECTION 3 MATERIAL 


The general approach common to panel methods was 
described in section 3. Here, a brief description is pre- 
sented of how that approach is actually implemented in 
the PAN AIR code. 

Recall that PAN AIR breaks each panel into eight tri- 
angular subpanels to maintain doublet-strength continuity 
within a network of panels. Consequently, separate linear 
source and quadratic doublet distributions are used over 
each subpanel, instead of there being a single source and 
single doublet distribution over the entire panel. The 
source distributions are chosen to be linear polynomials 

<j(^il) = ao+CT n Ti + a^ (35) 

and the doublet distributions are chosen to be quadratic 
polynomials 

ti ^ 

^.•n) = no+ ^ + ^mTF 

(36) 

where £ and t| are the in-plane coordinates of a local 
(C/n*C) coordinate system associated with each subpanel. 

Panel and Network Singularity Parameters 

The linear sources contribute three unknown coeffi- 
cients (Oo&fyGj]) per subpanel, and the quadratic doublets 
contribute six additional unknown coefficients per sub- 
panel, giving a total of 72 unknowns per panel. To solve a 
system of equations involving 72M unknowns, where M 
is the number of panels, would be extremely inefficient. A 
better approach is to relate the 72 subpanel coefficients to 
values of the source-doublet strengths at a small number 
of discrete points. PAN AIR uses two sets of points: one 
is used to define panel singularity parameters, and the 
other is used to define the network singularity parameters. 
The panel singularity parameters are used to evaluate 
panel influence-coefficients, and to relate the 72 coeffi- 
cients to the network singularity parameters {X} appear- 
ing in equation (17) in section 3: 

[AIC]{X} = {b} 


The eight subpanels of a panel are shown in fig- 
ure 34, along with the network singularity-parameter 
points (locations). Both source and doublet singularity 
parameters are defined at panel center points. Only 
doublet singularity parameters are defined along the net- 
work edges; these are used to enforce doublet strength 
continuity across network abutments. We will not discuss 
this aspect of PAN AIR; we will only discuss the 
relationships between the 72 coefficients and the panel 
and network singularity parameters. 

The panel singularity-parameter points are shown in 
the lower portion of figure 35. There are five points for 
the source and nine points for the doublet. The panel 
source and doublet singularity parameters are the values 
of the source and doublet strengths at these points; they 

are labeled o\ as, and pi,...,|i 9 , respectively. These 

panel source and doublet singularity parameters are 
expressed in terms of 9 network source and 21 doublet 
singularity parameters, {X s } and {X D }, respectively, 
located at the panel center point and at the neighboring 
panel center points indicated by the asterisks in the upper 
portion of figure 35. (When a panel is nearer a network 
edge than is shown in fig. 35, neighboring singularity 
parameters that are different from those shown are used.) 
This is done through what the PAN AIR theory document 
calls source and doublet “outer” spline matrices [B s ] and 
[B^], respectively (ref. 11). There is one each of these 
two matrices for each panel. 

To relate the panel singularity-strength parameters 
(shown in the lower half of fig. 35) to the three source and 
six doublet coefficients appearing in equations (35) and 
(36), PAN AIR uses source and doublet “subpanel” spline 
matrices [SPSPL s ] and [SPSPL D ], respectively. There is 
one each of these two matrices for each subpanel. 

The actual construction of the outer and subpanel 
spline matrices is described in appendix I of reference 11. 
Here, we simply show how they are used. For the sources, 
we can now use equation (35) and the matrices [SPSPL s ] 
and [B s ] to obtain, for subpanel k. 
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components at an arbitrary point P in terms of the panel 
singularity parameters (five for the source, nine for the 
doublet). That is, a relationship of the form 
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(39) 


is determined, where [PIC S ] and [PIC D ] are panel 
influence-coefficient matrices for the sources and dou- 
blets, respectively. This is done as follows. To determine 
the potential, we start with equation (9): 



+ p(Q)n- V q-^- 



(40) 


Similarly, equation (36) gives, for the doublet distribution 
on subpanel k, 

^•■n) k 


Ho 

H 


1 % T, 






Using equation (37), the potential owing to the panel 
source distributions is 

5 iI spsplS } 

(41) 
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where represents one of the panel’s eight subpanels. 
The subpanel spline matrix and the vector of panel source 
singularity parameters are not functions of the local sub- 
panel coordinates and can be extracted from the integral, 
giving 


* S <P)^X JJ ik 1 * nJd^dT! [SPSPL 5 ] 
k=1 L J A k J 
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CT 5. 

(42) 


(38) 

Panel Influence Coefficient Matrices 

Equations (37) and (38) relate the source and doublet 
distributions, respectively, of a subpanel, to the panel and 
network singularity parameters. Equations (9) and (1 1) are 
then used to obtain the perturbation potential and velocity 


Comparing equations (42) and (39), we see that the under- 
lined portion of equation (42) gives the elements of the 
first row of the [PIC S ] matrix. The same approach is used 
to obtain the other rows of both matrices appearing in 
equation (40). Of course all the surface integrals similar to 
the one appearing in equation (42) must still be evaluated. 
This step is described in appendix B. 
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Assembled Influence Coefficients 


(45b) 


Equation (39) gives the contribution of a single panel 
to the perturbation potential and velocity at a point. To 
arrive at the equation [AIC] {X.} = (b), we have to assem- 
ble the effects of each panel. This is simply done with 
superposition, which is allowed since we are solving dis- 
cretized equations corresponding to a linear partial differ- 
ential equation. 


r*(p)i 

’L<t>icj' 

jv(P)J 
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or 


lxN 


3xN 


<t>(P) = [<t>ICj{X} , v(P) = [VIC] {X} 


(45c) 


The first step is to use the outer spline matrices to 
express equation (39) in terms of the network singularity 
parameters, namely 


f«P)l 

v(P)j 


= [PIC S ][B S ]{X S } + [pIC d ][b d ][x d } (43) 


[ic s 


panel 


ic D 


panel 


The next step is to move to another panel, repeat the 
construction of equation (43), and superimpose the 
results. An example of this, for the source contribution to 
the potential at point P owing to two adjacent panels, is 
shown in figure 36. Applying this process to all the panels 
of each network then produces an equation of the form 


The latter form is what we need to use with the boundary 
conditions to finally arrive at [AIC] {X} = {b} . 

The Aerodynamic Influence Coefficient Matrix 
Equation 

Recall from section 4 that the velocity boundary con- 
dition for an impermeable panel is 

(v„ + v).n«0 (46) 


This is a specific example of the more general boundary 
condition, for point P, namely 


H 

lv(P)J 


4xN s 4xN d 
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L J I iD 


c«t>(P) + (an + t ) • v(P) = b(P) (47) 

(44) 

where c and a are arbitrary constants, and t is a vector 
tangent to the panel. 


where Ns = the number of network source singularity 
parameters, Nq = the number of network doublet singu- 
larity parameters, and N = Ns + Nd is the total number of 
singularity parameters for all the networks of the configu- 
ration (except for those involved with doublet matching 
across network abutments, which we have ignored). 
Equation (44) gives the effect of all panels on the 
(perturbation) potential and velocity at point P, in terms of 
all the network singularity parameters. Next, further parti- 
tion equation (44) as 



_4xN s 

4>IC S 
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(45a) 


Equation (46) is obtained from (47) by setting c = 0, 
a = 1 , t =0, and b(P) = -V„ ■ n . Equation (47) can be 
written as 


3xN 

c<J>(P ) + (an + t} T {v(P)} = b(P) (48) 


Then, using equations (45c), we obtain 
lxN 1x3 3xN 


c|_<t>ICj + {an + t} T [WIC]\{X] = b(P) (49) 


The underlined portion of equation (49) is a row of what 
we have been calling the [AIC] matrix, since if N bound- 
ary conditions are imposed, we will have N equations like 
(49) relating the N unknown singularity parameters; that 
is 
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[AIC]{A.} = {b} (50) 

Two issues have been sidestepped in the above dis- 
cussion. First is the question of doublet matching at net- 
work abutments. The reader is referred to the PAN AIR 
theory document (appendix F of ref. 1 1) for information 
on the doublet matching. Second is the question of how to 
analytically evaluate the surface integrals of the type 
appearing in equation (42). The complete story is con- 
tained in appendix J of reference 1 1 (approximately 
200 pages). An abbreviated version, restricted to incom- 
pressible flow, is given herein in appendix B. Finally, 


there is a third point which was mentioned in the Sum- 
mary of Methods portion of section 4. For supersonic 
flow, the outer doublet spline matrices have to be 
“forward-weighted” to obtain stable results at high Mach 
numbers. The weighting refers to the weights used in a 
least-squares construction of the outer spline matrix. The 
forward-weighting obtains more information from 
upstream points than from downstream points. This is 
somewhat analogous to the use of forward-differencing in 
supersonic zones by finite-difference codes. This forward- 
weighting scheme is described in section 1. 1.2.4 of 
reference 11. 
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APPENDIX B 


EVALUATION OF INFLUENCE-COEFFICIENT SURFACE INTEGRALS 


In appendix A, equation (40) contained the 
expression 

<K p ) = Y JJ s “TT^ dS Q < 51 > 


and, as shown in figure 37, X denotes the flat panel, 
Q(£,ti, 0) is a variable point of the panel, P(x,y,z) is the 
influenced point, and h is the height of P above or below 
the panel (h is positive above the panel, and negative 
below it). The x-component of velocity at P is 


for the potential at P owing to a source distribution over a 
panel. The equations following equation (40) were for the 
PAN AIR implementation wherein a panel is subdivided 
into eight triangular subpanels. Here, we focus on a flat 
quadrilateral (or triangular) panel wherein a single source 
distribution covers the entire panel. This eliminates the 
subpanel spline matrices of appendix A, thereby simplify- 
ing the discussion (but it also prevents the doublet 
strength from being continuous). We also restrict our dis- 
cussion to incompressible flow. The approach is that of 
sections D.5 and G.l in reference 33. The general scheme 
is to integrate the surface integrals by parts (where possi- 
ble); this reduces, by one, the singularity order of the 
resulting surface integrals, and introduces line integrals 
along the panel edges. The resulting expressions for the 
potential and velocity field are expressed below in terms 
of surface integrals denoted H(M,N,K), and line integrals 
denoted F(M,N,K), where M,N,K are integer exponents 
appearing in the integrands. We then show how H(l,l,3) 
and F(l,l,l) are evaluated analytically. These two inte- 
grals are the fundamental ones because, as shown in refer- 
ence 33, all the remaining integrals can be computed 
recursively in terms of these two. It is H(1 ,1,3) that causes 
the jump properties discussed earlier. The jump behavior 
will also be derived. 

Velocity Owing to a Flat, Linear Source Panel 

For a flat panel with a linear source distribution and 
incompressible flow, equation (51) becomes 


= = d ^ d11 " 4 ^ r ) CTd ^ dT1 


(54) 


Noting that 


3(1 /R) 3(1 /R) 

5x 3£ 


equation (54) is written as 

<56) 


This form allows us to integrate by parts, using the 
formula shown in figure 38 (this result can be derived 
from the equations on p. 499 of ref. 63), to obtain 



xi 


.R-OV^- 



1 3cr JE . 


(57a) 


where the sum is over the four panel edges. Similarly, the 
y-component of velocity is 



(57b) 



R(£,Ti;x,y,z) 


d^dri 


where 




and 


(52) 


(53) 


The z-component of velocity is 

3 (1 /R) 
3z 


= 3 £ = __Lff 
1 3z 4 jiJJ l 


CTd^dq 


ad£dq 


(57c) 


R= (S -*) 2 + (il-y) 2 + z 2 ' =[r 2 + h 2 ] 1/2 


We now collect terms: 
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A A A 

v= v x i + v y j + v z k 

* i + 'R’^ j+0 k ] <i5 *' 

+ iJJ £ [»- 5 + » 

and for convenience, write the above as the vector 
equation 


This is done to emphasize that these terms are indepen- 
dent of % and q. Also define the general line integral 
F(M,N,K) as 


F(M, 


,N,K) = J^- 


(^-x) M Vj-y) 


N-l 


,K 


-d^ 


(62) 


Then, equation (60) becomes 

(vb.Vjj.O) r 

5 b = 1 {o(x.y)F(l, 1, 1) + a x F(2, 1, 1) 


+a, 


r F(l,2,l)J 


(63) 


v = v A + v B + v c (59) 

where the three terms on the right-hand side of equation 
(59) are simply shorthand notations for the three lines, 
respectively, on the right-hand side of equation (58). 

The next step is to express the various integrals 
appearing in equation (58) in terms of (£ - x) and (q - y). 
This is done to put the integrals in a form that can be inte- 
grated analytically. To demonstrate how this is done, we 
will show the process for v B and v c . The process is the 
same for v A , 


At this point let us pause and reflect on what is 
known and unknown. In principle, for an arbitrary point P 
the line integrals F(M,N,K) can be evaluated for each 
edge of a panel and are, therefore, known. The coeffi- 
cients a 0 , (T£, and 0^, are unknown. They are, however, 
related to the source singularity parameters through an 
outer spline matrix as discussed in appendix A. (The 
actual spline matrix used in ref. 33 is given by eq. (B-7) of 
that reference.) Thus, once the singularity parameters (X) 
are determined from [AIC] {X} = (b) , everything in equa- 
tion (63) is known, and the velocity at P (from the v B 
term) can be computed. 


Starting with v B , the line integral contribution from 
edge / is 


The procedure for the surface integrals is similar, as 
we illustrate for the contribution of v c to equation (58): 




. K- v n-°)fo J( 

V B = -^T-J / R d/ 

Km) f / , \ i 

=1 ^ir^J / K +a ^ + MR 

= jj( q ° + + <v) + - *) 

+tf n ( T i - y)]i' d/ 


(60) 


Now define 


a(x,y)sa 0 +o^x + o n y 

a x = a 4 


(61) 


CTy-a^ 



(64) 

where equations (61) have been used. Now defining the 
general surface integrals H(M,N,K) as 

H(M,N,K) = jj r k^ ~ ~ <^ dl l ( 65 > 

permits us to write equation (64) as 
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VC = {a(x, y) • hH(l, 1. 3) 

+a x • hH(2,l,3)+a y ■ hH(l,2,3)J (66) 

If we had analytic expressions for the H(M,N,K) integrals, 
the v c contribution to the velocity field could be com- 
puted (in terms of the three unknown source coefficients) 
from equation (66). A similar expression is likewise 
obtained for the v A contribution. 


E2 = £(0,0-l)[(n x v^ + 4 y v n )F(l, 1, 1) 

1 

+(4xxV 4 + n xy v n )F(2,l,l) 

+ (^*y v $ + %y v n) F ( 1 - 2 > 1 )] (68b) 

and 


Velocity Owing to a Flat, Quadratic Doublet Panel 


n(S.'n) = v-o + + 


Before showing how the fundamental integrals 
H(l,l,3) and F(l,l,l) are evaluated, let us first discuss 
some of the properties of the various integrals (in practice, 
one only knows these properties after evaluating and 
studying the integrals). A more complete discussion 
appears in appendix D of reference 33. We will do this by 
first giving the complete equation for the velocity owing 
to a flat panel with the quadratic doublet distribution 
given by equation (36), and then describing the signifi- 
cance of the various terms. This equation is obtained from 
equations (D.130) through (D.140) of reference 33; the 
derivation follows the procedures described above for the 
source panel: 

4itv= (0,0,1) •[ji xx + n yy ]-H(l, 1,1) 

+[(m,H y ,o)- hH(l,l,3) • (Hxx» o) • hH(2, 1, 3) 

+(m y .Hyy.o)h ■ H(l,2,3)] + El + E2 


(67) 


where 


+ 7*Vi t 1 2 ’ M) g2; 

n(x. y) = Vo + + n n y + \ 


1 2 

+ I ^y 


(69) 


Mx(x.y) = ^ + ^x + ^, 1 y 
Hx* = 
and so forth. 

The first two terms in equation (67) arise from the 
surface integrals produced by the integration by parts — 
the last two terms arise from the line integrals. The first 
and last terms give velocity contributions that are in a 
direction normal to the panel; the second gives contribu- 
tions that are in directions parallel to the panel. The third 
term is the "line-vortex” contribution and produces a 
velocity that is in the direction of / x g, which is shown in 
figure 39. 


hv^ -h\ n -a 
g ’ g ’T 

unit vector in 
direction of /®g 

+H X F(2,1,3) + |i y F(l, 2,3) + £u„F(3, 1.3) 




+H xy F(2,2,3) + ^ yy F(l,l,3) 


(68a) 


Both F(l,l,l) and F(l,l,3) would be singular if they 
were to be evaluated at a point of the panel edge. In prac- 
tice, this problem is avoided by having edge control points 
(if used at all) located slightly interior to the panel. For 
points near a panel edge, F(l,l,l) is the weaker singular- 
ity, varying as ln(l/g*), whereas F(l,l,3) varies as 1/g 
(like a line vortex). 

The quantity hH(l,l,3) in the second term is what 
causes a panel’s jump properties. Here, we see that it 
multiplies the in-plane first derivatives of the doublet 
strength; this is what causes the tangential velocities of a 
doublet panel to be discontinuous across the panel. Note 
that this property would not be predicted if the panel dou- 
blet strength were assumed to be constant. 
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The surface integral H(l,l,3) is singular if computed 
at an interior point of the panel, but hH(l,l,3) is well 
behaved (just having the jump behavior). For a point on 
the panel edge, hH(l,l,3) depends on the angle at which 
the point is approached. This nonuniqueness is avoided as 
done for the F integral singularity behavior, that is, by 
using edge control points withdrawn slightly into the inte- 
rior of the panel. The behavior of the doublet velocity 
near a panel edge can be used to impose the Kutta condi- 
tion (see fig. 10 in ref. 9). 

If the doublet strength is continuous from panel to 
panel (as in the PAN AIR, HISSS, MCAERO, and 
QUADPAN codes), the line-vortex contributions from 
adjacent panel edges will be equal and opposite and will 
cancel, and thus will not produce a resultant velocity field. 
(In general, a correct potential flow solution to flow about 
a general shape does not contain concentrated line- vortex- 
induced velocities.) 

Ironically, if a constant-doublet-strength panel is 
used, the linear and quadratic terms in equation (67) dis- 
appear, and the only nonvanishing term is the spurious 
line-vortex term associated with F(l,l,3); that is 

4 

* = ^ll/® g40gF(1,u) (70) 

I 

If the above expression is evaluated for a specific 
panel geometry, the velocity field obtained will be the 
same as that produced by four line-vortex filaments all of 
strength T = Po, placed head-to-tail around the panel 
perimeter. That is, a constant-strength-doublet panel is 
equivalent to a ring vortex. Although many subsonic-only 
codes use the line-vortex model (and often get away with 
it numerically), such an approach causes numerical disas- 
ter in supersonic flows. 

Evaluation of H(l,l>3) 

The surface integrals H(M,N,K) were defined by 
equation (65). The fundamental integral H( 1,1,3) is then 

h(u ' 3) -1L^ <7i> 

This integral can be evaluated by making two changes of 
variable. The first is to use polar coordinates in the plane 
of the panel, as indicated in figure 40, giving 



The upper limit r extends to the boundary of Z; h is a con- 
stant as far as the integration is concerned; and i is the 
panel comer point number. 


Next, the integration on <|> is changed to a line integral 
along the panel edges. Referring to figures 40 and 41, 
equation (72) becomes 



where 


g 2 = a 2 + h 2 


(73b) 


and where a is a different value for each edge L. The 
geometrical significance of g is shown in Figure 42. The 
indefinite integrals of the two terms appearing in equa- 
tion (73) can be found in tables (e.g., on p. 49, line 3 of 
ref. 64 for the second term, called I2 below). The result is 



where 


T f /i+1 d t 1 , -if t \ *i+i 

''“J,, FTF' 1 ® N<, 


and 


. . r' 1 *' & ±_ 


(/ 2 + 5 2 )/<W 


1 ...-1 


( 


TaTThi 


tan 


I hi/ 


lal^ 2 + g 2 


*i + i 

A 


(74) 


(75a) 


(75b) 


Thus, we can write 
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( 76 ) 


The Jump Property of hH(l,l>3) 


which is equation (21), except total rather than perturba- 
tion velocity appears in that equation. Both equations are 
correct, since the total and perturbation velocities differ 
only by the free- stream value, which is a constant. 

The Line Integral F(l,l,l) 


From equation (62), and figure 41, the fundamental 
line integral F( 1,1,1), for edge L, is 



This can be put in several forms, for example. 



Analytically, these three forms are identical. How- 
ever, depending on the location of point P, the different 
forms produce different numerical results with a com- 
puter’s finite precision arithmetic. Referring to figure 41, 
the following choice is recommended in reference 33 as 
the best for numerical accuracy. 

When Use 


o F[(1,U) 
/,,/ 2 <0 F 2 (1,U) 

*!<(U 2 >0 Fj(l,l,l) 
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APPENDIX C 


IMPROPER INTEGRALS AND THEIR FINITE PART 


The velocity equations for the source and doublet 
distributions both contain the term hH(l,l,3) (see 
eqs. (66) and (67), respectively). The integral contained 
therein was evaluated for an arbitrary point P at (x,y,z). In 
particular, z was nonzero. For P being an interior point of 
the panel, z = h = 0, we found that hH(l,l,3) is generally 
well behaved. It has different values on opposite sides of 
the panel, but it is not singular. 

If z = 0 had been used before evaluating the integral, 
the resulting integrand would have been simpler, but it 
would also have been singular. It turns out that singular 
(improper) integrals obtained in this manner can actually 
be correctly evaluated. This behavior is illustrated with 
the following two-dimensional example. 

We are going to compute the two-dimensional veloc- 
ity field owing to a doublet sheet (panel) in the x-y plane. 
The panel extends to infinity in the positive and negative 
y-directions, the doublet strength varies only in the 
x-direction, and the doublet axis is in the z-direction 
(normal to the panel). 


We start with the expression for the potential at 
P(x,z) owing to a line doublet extending from y = \o 
y = +oo, passing through the point x = z = 0. The strength 
per unit length of the y-direction is p, and the doublet 
axis is aligned with the z-direction: 



(89) 


If instead of a single line doublet, we have an 
x-distribution of them located between x = a and x = b, of 
strength p,(xi) per unit length of the x-direction, the poten- 
tial at (x,z) is 




Xi=a(x-X 1 ) +Z‘ 


-dx\ 


(92) 


We are going to evaluate the expression for w(x,0) in two 
different ways. The first is to do the integration first, and 
then set z = 0 to find w(x,0). This is how we did it in 
appendix B for hH(l,l,3). The second way is to do the 
integration after setting z = 0; this will produce an 
improper integral containing a nonintegrable singularity. 

Integration Before Setting z = 0 

Performing the differentiation in equation (92) gives 




(93) 


If the velocity is wanted at a point x in the interval 
a < x < b, the term (x - xi) will become zero as the 
dummy variable of integration, xi, passes through the 
point x. However, the integrand will remain finite if we do 
not let z = 0 at this stage. 

Performing the integration yields 

Xl=b 


Xi=a 


w(x,z) = ^l 


X-Xj 


(x-Xj) + TT 




K*i) z 


x 1 =a(x-X 1 ) 2 + Z 2 




(90) 


Letting p.(xi) = (io, a constant, we have for the potential 
and z-component of velocity, the expressions 

<9,) 


Now we let z = 0, obtaining 



-1 1 ' 
x-a x - b 


(95) 


This expression gives the z-component of velocity at any 
point along the x-axis; it is the same as that owing to a 
pair of concentrated line vortices (each of which produces 
a velocity v = r/( 27 tr)) located at x = a and at x = b, whose 
strengths are equal and opposite, and of magnitude 
T= Ho. This is a confirmation of our earlier conclusion 
that the velocity field owing to a constant-strength doublet 
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panel is the same as that produced by line-vortex fila- 
ments (having the same strength as the doublet) being 
placed along the panel edges. 

Integration After Setting z = 0 


evaluate the nonsingular integral (the one containing the 
nonzero z), but it may be possible to find the indefinite 
integral when z = 0. Also, even if the nonsingular problem 
can be evaluated, it is usually simpler to evaluate the sin- 
gular problem, as it was in our example. 


Here, we set z = 0 in equation (93), yielding the 
integral 



(%) 


This is an improper integral, since if x is wanted at a point 
between a and b, the denominator is zero at xi = x. 
Because the integrand becomes infinite, this integral has 
no meaning in the ordinary sense. However, let us ignore 
this problem, and proceed as if the integrand were not 
singular. We then obtain 


i i*i = b _^ 


w ( x *°) 2*x-Xil Xl _ lx 


-1 __ 
x-a x 


1 ' 
-b 


(97) 


Amazingly, we have obtained the correct answer, that 
is, we have obtained equation (95) again, even though the 
integrand is singular. This result is not a fluke. Refer- 
ence 65 (p. 1 1) shows that, if 


There may be cases in which even the indefinite inte- 
gral for the singular integral is not known. In such cases 
the finite part can be obtained numerically. One way is as 
follows. If 


r b f( Xl )d Xl 

'*!=« (*- X l) 2 


( 100 ) 


is the improper integral obtained by setting z = 0 before 
doing the integration, the finite part can be computed from 
the formula 

r b fylfot 

Jx!=. (X-Xi ) 2 

.uJpibfe+r'’ ib&.m 

e- + 0 [Jx 1 =a(x-X]) Jxi=x+e(x-Xj) 


( 101 ) 


*> 

J (x-x,) 

is the indefinite integral, for example, G = l/(x - xi) in 
our simple example, then the correct finite value for the 
improper integral (that results from setting z = 0 before 
performing the integration) is 

f b 7 ^%r = G(x,b)-G(x,a) n = 0,1,2,... (99) 

^.(x-Xj) 

The expression on the left-hand side is merely notation 
used to indicate that the integrand is singular, but that the 
integral does in fact have a unique finite value. This value 
is usually called the "principal value” or the "finite part” 
of what appears to be an integral having an infinite value. 
Equation (99) says that this finite value can be computed 
from the indefinite integral by evaluating it at the two 
end-points a and b, as we did in our example. If the indef- 
inite integral contains logarithmic terms, the absolute 
value of the arguments must be taken. 

The reason for dealing with the finite part concept is 
that unlike our simple example, it may not be possible to 


where e is a small number. For a fixed value of x and e, 
the two integrals on the right-hand side of equation (101) 
can be evaluated numerically. Their sum will be a "large” 
number, from which is subtracted another large number 
given by the last term. The limiting value of this differ- 
ence, as e goes to zero, exists and is the finite part. Equa- 
tion (101) is a special case of the more general formula 
given by equation (13) of reference 64. Another numerical 
approach is given in reference 66 (pp. 42-47). 

Contrary to what one might assume based on the 
above discussion, performing the integration before set- 
ting z = 0 does not always eliminate the appearance of 
improper integrals. A counter-example appears in refer- 
ence 33 (appendix D.2), where a curved panel is exam- 
ined. In that case, the integration-by-parts step we went 
through in appendix B is not done. The resulting integrals 
are all expressed in terms of a nonzero z, but as z 
approaches zero some of the integrals become improper. 
However, the collection of all the improper integrals sums 
identically to zero, and the improper integrals can conse- 
quently be ignored. 

Additional information on improper integrals can be 
found in reference 37 (pp. 15-18) and in reference 11 
(appendix J). 
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TABLE 1.- COMPARISON OF FORMULATIONS 



Approximate 


Singularity type j 


Code name 

date 

introduced 

Panel 

geometry 

a 

n 

Y 

Reference 


Moo, subsonic only 


Lower-Order 

Hess 

1965 

Flat 

Constant 


Linear 

2,38 

MCAERO 

1980 

Piecewise 

flat 0 

Constant 

Quadratic 42 


17 

SOUSSA 

1980 

Hyperbolic 

paraboloid^ 

Constant 

Constant 


8 

VS AERO 

1982 

Flat 

Constant 

Constant 


18,19 

LEV 

1979 

Hyperbolic 

paraboloid 0 

Linear 

Quadratic 42 


27 

QUADPAN 

1983 

Flat 

Constant 

Constant 


20,21 

Higher-Order 

Hess 

1985 

Parabolic 

Linear 


Linear 

32 


Moo, subsonic/supersonic 


Woodward- 

Carmichael 

1966 

Flat 

Line distributions 
for bodies 

Constant 

4,5 

USSAERO 

1973 

Flat 

Constant 


Linear in x 
Constant in y 

7 

PAN AIR 

1981 

Piecewise 

flal° 

Linear 

Quadratic 0 


9-15, 

36,37 

HISSS 

1984 

Piecewise 

flat 0 

Linear 

Quadratic 0 


16,42 

QUADPAN 

1985 

Piecewise 

flat 0 

Constant 

Linear 42 


b 


''Continuous panel edge geometry or doublet strength. 

J’From October 1985 AIAA Workshop on Aerodynamic Analysis Using Panel Methods. 


35 






KUTTA CONDITION 



WAKE 


a) SUBSONIC (1-M «£) 0 XX + <£ yy + tf> zz = 0 



b) SUPERSONIC ) </> xx - 0 yy - 4> zz - 0 

Figure 3 - Panel methods predict linear potential flow by solving the Prandtl-Glauert equation. The weak shock solutions 
(small 0) of shock-expansion theory are predicted and yield wave drag. 


MODIFIED BODY 



Figure 4.- Two methods of boundary-layer simulation (not to scale) (from reference 25). (a) Surface displacement; 
panels at modified body, (b) Surface blowing; panels at actual body. 
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CENTER OF PANEL i 



KNOWN SOURCE AND DOUBLET SOLUTIONS 
ARE FUNCTIONS OF R AND 



7 0j =// [a(xj) x K 0 +m(x|)xK^] dSj 

/ 5 / ~T 

/ UNKNOWN SOURCE AND DOUBLET STRENGTH 
/ DISTRIBUTIONS AT PANEL j 

VELOCITY POTENTIAL AT POINT i,^| = V 0j IS THE VELOCITY AT POINT i 


Figure 7 - The flow at center of panel i due to source and doublet distributions at panel j. 


HORIZONTAL AND 
VERTICAL TAIL 

MIDBODY AND NETWORKS 



Figure 8 - Networks of surface panels. Each panel is a surface distribution of sources and doublets whose strengths 
adjusted so as to satisfy the boundary conditions; for example, to make the flow tangent to the surfaces. 



tg differei 
for visual 






n 


(a) 

Figure 10.- M x 1 



CORNER POINTS 1 AND 3 ARE BELOW THE 
PLANE OF THE FLAT CENTRAL PORTION AND 

CORNER POINTS 2 AND 4 ARE ABOVE THE 
PLANE OF THE FLAT CENTRAL PORTION 


STRAIGHT LINE PANEL EDGES 




(a) (b) 



Figure 12.- Continuous panels constructed from four arbitrary points, (a) Piecewise flat; QUADPAN. (b) Piecewise flat; 
PAN AIR, HISSS, MCAERO. (c) Hyperbolic paraboloid; SOUSSA, LEV. (d) Figure 12(c) corresponds to the 
twisted shape a flat flexible plate takes when loaded as shown. 
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• TWO SIDES 
OF A PANEL 



• SOURCE 



• DOUBLET 



Figure 13.- The jump properties of source and doublet panels; incompressible flow. 


ACTUAL 



SURFACE 



Figure 14 - The spurious line vortex that results from having a discontinuous doublet-strength distribution. 
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THEORY 

O EXPERIMENT 




Figure 15- Effect of discontinuous doublet strength in supersonic flow: Moo = 2.05, a = 2°, 3%-thick biconvex airfoil. 


SLOPE OF SOLID SURFACES 
MUST BE SMALLER THAN THAT 
OF THE MACH CONE 



Figure 16 - Supersonic flow modeling constraints. 
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Figure 17.- Two-dimensional pressure distribution for supersonic flow over a wedge. 



Figure 18 - Jump in mass-flux and velocity vectors across Mach line from wedge apex. 
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\ 



Figure 19 - Cone pressure solutions to Euler and Prandtl-Glauert equations at supersonic Mach numbers. 



Figure 20 - Alternative ways of representing a solid surface with source/doublet panels, (a) Direct approach: <(> int = 0, 

W ext • n = 0 . (b) Indirect approach: 4>int = 0,<r = -V^ • n . 
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Figure 21.- PAN AIR model of supersonic Fighter concept. 
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O WIND TUNNEL DATA 
O PANAIR ISENTROPIC RULE 
□ PANAIR SECOND-ORDER RULE 
A USSAERO-B 


Figure 23 - Lift and moment coefficients for configuration of figure 21: M»o = 2.2. (a) Lift coefficient, (b) Pitching 

moment coefficient 








' ^ 7 ^ 


-X 



M 


TT*~ 

yl 



Figure 24.- A wing-canard VSTOL concept. 



dimensionless spanwise 

CIRCULATION DISTRIBUTION 



UNALIGNED WITH 



ALIGNED WITH 


Figure 25 - Canard-wake models: a = 10°, root section. 



WING SPANWISE POSITION, rj- 


Figure 26.- Effect of canard- wake position on spanwise circulation distribution. 


50 




Figure 27- Paneling of Boeing 737-300 flaps-15 configuration. 



Figure 28 - Flaps-15 modeling, (a) Flap geometry, (b) Wakes. 
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• COMPUTATIONAL 
BUILDUP AT 1.2 V $ 



• NACELLE • WING • FULL AIRPLANE 

• HORIZONTAL @1.2V $ @1.2V S 

• VERTICAL 

• BODY 

METHOD — ►HANDBOOK 2D MULTIELEMENT PAN AIR 

VISCOUS 

(a) (b) 

Figure 29.- Predicted and measured drag at 1.2V S ; V s = stall speed, Cop = profile drag, = induced drag, FI = flaps-! 
setting, F15 = flaps- 15 setting, (a) Drag breakdown at two flap settings, (b) Aircraft lift-to-drag ratio. 
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Figure 33.- Detail of leading-edge paneling embedded within the rectangular grid. 


O LOCATION OF DOUBLET SINGULARITY 
PARAMETERS \P 

* LOCATION OF SOURCE AND DOUBLET 
SINGULARITY PARAMETERS Af , A 1 ? 

A? AND A? ARE THE UNKNOWNS, [AIC] { A ) = { b } 


A PANEL AND 




(a) (b) 

Figure 34.- A nine-panel network and its singularity parameter locations, (a) Network, panel, subpanels, (b) Location of 

network singularity parameters and corresponding control points. 
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SOURCES 


DOUBLETS 
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* 

* 
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1 'T' ' 

/ \ 

\i / 

* 



* 

* 

* 








* = THE 9 NETWORK SINGULARITY 
PARAMETERS A? ASSOCIATED WITH THE 
5 PANEL PARAMETERS 



* 

* 

* 


* 

* 

* 

* 

* 

* 

* « 

\ i / 

, v i 

• * 

* 

* 

* 

* 

* 

* 


* 

* 

* 



* = THE 21 NETWORK SINGULARITY 
PARAMETERS A? ASSOCIATED WITH THE 
9 PANEL PARAMETERS 



• =THE 5 PANEL-SINGULARITY 
PARAMETERS 


"5 



•=THE 9 PANEL-SINGULARITY 
PARAMETERS 


Figure 35.- Panel-singularity parameters and the associated network-singularity parameters. 
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12 x 1 




Figure 36 - The potential at P owing to the source distributions of panels 6 and 7, in terms of the associated network 

singularity parameters. 



Figure 37.- Flat panel £ with linear source <tq and influenced point P. 
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//f.jgdS = ^fgi'jdC - //fg.jdS 

i = 1 -*• v i = ^ 

' = 2 - v 2 = v n 




4 

//f -gd^dTj = Z / fg^jdC: - //fg.jd£d7? 
Z i =1 Cj 2 

j = EDGE NUMBER 


(b) 


Figure 38 - Two-dimensional integration by parts, (a) General form, (b) Form for a quadrilateral panel. 



Figure 39.- The vector / x g . See figure 41 for definition of a . 
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V'V 



ON PLANE OF I 

Figure 40.- Polar coordinates (r,4>) and edge coordinated^. 



FOR P TO THE 
LEFT OF L 


ALONG L WE HAVE 
r 2 = i 2 + e 2 


COS0 = 


vA 2 + a 2 


sgn< 6 — Lsiqn(a) 

\A 2 + a 2 

ta n<p = C/a 

d(t£n£ = sec 2 0 d0 = l 

dC dC a 


d0 = 


adC 


e 2 + a - 2 


Figure 41- Relationships between polar coordinates (r,<(») and the edge coordinate l . 
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